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NOMENCLATURE 


A 

C p specific heat capacity (Btu/lb-°F) 

E dependent variable, H or T 

H enthalpy (Btu/lb) 

h heat transfer coefficient (Btu/ft 2 -hr-°F) 

K thermal conductivity (Btu/ft-hr-°F) 

KW average thermal conductivity (Btu/f t-hr-°F) 

L the total length of the composite slab in the y-direction (ft) 
L f latent heat of fusion for ice (Btu/lb) 

3 

q rate of heat generation per unit volume (Btu/hr-ft ) 

q M rate of heat generation per unit area (Btu/hr-ft or Watts/in ) 

T temperature (°F) 

t time variable (hr) 

x space coordinate in x-direction (ft) 

x dimensionless space coordinate in x-direction 

y space coordinate in y-direction (ft) 

dimensionless space coordinate in y-direction 


Greek Letters; 


a 

At 

Ax 


thermal diffusivity (ft /hr) 


Iii 


size of time step (hr) 
dimensionless grid spacing 


NOMENCLATURE (CONTINUED) 


Greek Letters; 

Ay dimensionless grid spacing 

0 dimensionless temperature 

p density (lb/ft^) 

Subscripts: 

b^ inner ambient boundary 

b 2 outer ambient boundary 

g the position in the center line of the gap 

h the position in the center line of the heater 

1 grid point in x-direction 

j layer in composite blade 

k grid point in y-direction 

1 liquid (water) 

lm liquid at the melting point 

m melting point 

s solid (ice) 

sm solid at the melting point 

Superscripts: 

A evaluated at the previous time level 

o evaluated halfway between the previous and present time levels 

N value from previous iteration 

N+l value from current iteration 
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I. INTRODUCTION 


Many practical transient heat conduction problems involve phase 
changes, e.g., the melting of a solid or the solidification of a liquid. 
These types of phase change problems are known as Stefan problems, after 
the early work done by Stefan (1889), and they have recently received 
considerable attention. Two of the oldest examples of this type of 
problem are the formation of ice on a solid body and the solidification 
of lava streams. Other applications include ablation of missile skins 
under aerodynamic heating, casting of metals and glass in molds, freeze 
drying of foodstuffs, heat exchangers using cryogenics, and the freezing 
and thawing of soil. The number and diversity of melting/ freezing 
applications has motivated much research, which has been extensively 
reported in the literature [1]. 

The formation of ice on aircraft components has a considerable 
effect on aircraft performance as it increases drag and decreases lift. 
Thus, an aircraft must be designed with the necessary equipment required 
for ice removal or prevention. Basically, ice protection systems can 
be classified as either anti-icing or de-icing. 

The anti-icing (or ice prevention) principle involves the prevention 
of ice formation on the protected area at all times. Typical methods 
of anti-icing employ the use of chemical systems (freezing point depres- 
sants) or the passage of hot, compressed bleed air from the engine 


1 


through passages below the surface on which ice formation is to be 
prevented. Anti-icing methods usually utilize large amounts of energy. 
On the other hand, the de-icing (or ice removal) principle involves Che 
periodic removal of a layer of ice by mechanical or thermal means. This 
is accomplished by destroying the adhesion strength between the ice and 
the aircraft surfaces. Various de-icing techniques that have been 
Investigated and reported upon [2, 3) are pneumatic boots and thermal 
techniques, the latter comprising cyclic heating of discrete elements 
by hot fluid or electrothermal means. For de-icing systems, the energy 
requirements are significantly less than those required for anti-icing 
systems. 

The choice between anti- icing and de-icing techniques depends upon 
the component geometry and its function, and upon other factors such as 
the energy source utilized, the sensitivity of the component and the 
total energy available for ice protection. The advantages and pitfalls 
of anti-icing and de-icing systems have been reviewed by Baliga [4]. 
Stallabrass [5) has concluded from his experimental studies that the 
electrothermal method has the most advantages as a de-icing mechanism, 
although it does have maintainability/ reliability problems. Werner [3] 
also reported that the electrothermal de-icing technique is the method 
most commonly used, and it has been applied to both fixed and rotary 
wing aircraft. For these reasons, the electrothermal means of de-icing 
will be the major topic of the present work. 

The electrothermal de-icing principle was first successfully , 

applied to propeller blades. The purpose of this thermal technique is 
to raise the composite blade surface-ice interface temperature to above 
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the melting point of Ice, resulting In a very thin liquid interface 
layer which reduces the Ice adhesion to the blade surface. Aerodynamic 
and/or centrifugal forces can then sweep the unmelted ice from the 
surface. The basic two-dimensional electrothermal de-icing pad employed 
in this study is shown in Figure 1. It is essentially a composite body 
consisting of: (1) a metal substrate (aircraft blade)'; (2) an inner 

insulation layer; (3) a heater element; (A) an outer Insulation layer; 
and (3) an abrasion shield. A layer of ice is shown attached to the 
abrasion shield. 

The design of the heating elements is either spanwise on the blade 
to cause chordwise shedding, or chordwise to cause spanwise shedding. 

The heater element usually consists either of a woven mat of wires and 
glass fibers or of multiple strips of resistance ribbon. As pointed out 
by Werner (3], the key factors in the design of an electrical heating 
element for a cyclic de-icing system are: (1) the uniformity of the 

obtained surface temperature at the abrasion shield-ice interface; (2) 
the amount of Insulation required between the heater element and the ice 
surface; and (3) the amount of Insulation required between the heater 
element and the aircraft blade. Stallabrass [3] has mentioned that 
heater element configurations may vary considerably. One such configura- 
tion consists of an array of parallel wires, woven or knitted with nylon 
or glass filaments to a surface or plane heater mat; the wires might be 
0.02" in diameter, and spaced 24 to the inch, whereas metal ribbons 
usually have a thickness between 0.001" and 0.005". The gaps which exist 
between these heating elements can reduce the heater pad de-icing 
effectiveness, causing non-uniform melting of the ice. The gap width is 
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approximately 0.040" for metal ribbons and 0.080" for wires. 

The two Insulation layers, which usually consist of resin Impreg- 
nated glass cloth, serve to provide electrical insulation for the 
heating element. In order to direct more heat flow toward the Ice 
layer, it Is necessary to use a much greater thickness for the Inner 
insulation than for the outer insulation. A minimum thickness ratio of 
at least 2:1 has been suggested by Stallabrass f 5 J . The thickness of 
the outer insulation typically ranges from 0.010" to 0.020". 

The purposes of the abrasion shield, which is generally made of 
stainless steel, are: (1) to project the de-icer pad from rain erosion 

or sand and stone abrasion (especially when flying at high speeds); and 
(2) to diffuse the heat from the heater element so as to provide more 
uniform heating, thus minimizing cold spots above the heater gaps. The 
range of abrasion shield thickness is from 0.010" to 0.020". The material 
of construction of the substrate may vary widely depending upon the 
application. An aluminum alloy is considered in this study. 

A two-dimensional model of the de-icing pad is the subject of this 
investigation. The model, which represents a cross-section of the heater 
pad normal to the run of the heater elements, is shown in Figure 1. It 
can be seen in the figure that natural boundary lines have been formed 
by lines of symmetry of the centerlines of the heater strip and the gap 
between the heater strips. It is assumed that there is perfect adhesion 
between each layer and therefore no contact resistance to heat flow 
shall be introduced into the analysis. A point heat source or a finite 
thickness heater providing either a constant or a time dependent heat 
supply will be included in this study. A wide range of parametric 


4 


I 


studies have been conducted using the model to investigate such 
phenomena as the effect due to different heating functions on the rate 
and depth of ice melting, the effect of the abrasion shield toward 
evening out the temperature field due to non-uniform heating, and the 
effect of the gap width on the heater pad de-icing performance. A 
finite difference numerical analysis using the Enthalpy method to simu- 
late the phase change in the ice layer is employed in the model 
solution. 
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II. LITERATURE REVIEW 


A comprehensive review of the recent literature reveals that the 
study of electrothermal de-icer pads has either been ignored or that 
relevant information has not been published in the open literature. Of 
the investigations found, only the works of Stallabrass [5], Baliga [4] 
and Marano [6] have considered the effect of the phase change in the ice 
layer. All of the models have been one -dimensional except for that of 
Stallabrass, who proposed both one-dimensional and two-dimensional models. 
A two-dimensional electrothermal de-icer pad accounting for the phase 
change in the ice layer has been considered in the present study. 

Temperature distributions in composite bodies consisting of several 
layers of material have numerous applications in heat transfer problems, 
such as in rocket thrust chamber liners, fuel elements for nuclear 
reactors, and in reentry spacecraft. Several analytical methods may 

be used to obtain solutions for a limited number of layers with various 
boundary conditions. The Laplace transform technique, as presented in 
Carslaw and Jaeger [7], is the most common analytical technique used to 
solve transient heat conduction problems in composite bodies. In 
principle, this technique will yield the temperature distribution in any 
composite slab, but, in practice, if the number of layers'' exceeds two, 

i 

the inverse Laplace transform becomes quite difficult. Scott [ 8 J applied 
this method to determine the temperature distribution in a symmetrical 


five-component slab problem, and obtained a complicated solution which 
Is limited to constant physical properties within each slab. Goodman |9 J 
Introduced thiv adjoint -solution technique, which is a method of obtaining 
the solution to a large class of heat conduction problems in composite 
bodies from the solution of but one adjoint problem. The adjoint solu- 
tion arises from consideration of an auxiliary function with the 
application of Green's theorem to transform the problem from a volume to 
a surface integral. The primary disadvantages of this method are that 
the solution provides Just the interfacial temperatures and not the 
temperature profiles within each layer In the composite body, and that 
only steady state cases were evaluated. The extension of this method to 
handle the transient case has been done by Bouchillon [10]. Campbell 
[11J proposed a method analogous to electric transmission line theory to 
calculate the temperature profile in a composite body. Stallabrass [5] 
applied Campbell's method to check the accuracy of his numerical method. 

Tittle [12, 13] has shown that a one-dimensional multilayer problem 
can be solved by applying the Orthogonal-Expansion technique* This 
technique basically Involves constructing orthogonal sets from the solu- 
tions of each of the layers. An orthogonality factor, called the 
discontinuous weighting function, is used such that the resulting ortho- 
gonal set is applicable to the entire composite media. 

The disadvantage of using an analytical technique is that an 
excessive number of calculations must be done for each temperature desired. 
Hence, as the number of layers within the body increases, the calcula- 
tions become prohibitively cumbersome. In addition, analytical methods 
can not be applied when a phase change is considered. 
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The most commonly used numerical technique £or solving partial 
differential equations is the finite difference method, or, alternatively, 
the grid method. From a study of the recent literature, it appears that 
almost all of the numerical models for de-icer pads have used finU« 
difference methods. The method is considered approximate In that the 
continuous apace domain is replaced by a grid composed of nodal poi.its, 
and the derivative at a given point is represented by a derivative taken 
over a finite interval across the point. The differential equation is 
replaced by a system of difference equations. The difference equations 
are then solved by algebraic or arithmetic procedures. Using this 
method, the temperature at any or all of the nodal points within the 
composite body may be obtained at any or all time steps. The disadvan- 
tages of the analytical methods are thus overcome by the use of finite 
difference numerical methods. 

Two types of difference equations have previously been studied [14]: 

(1) explicit difference equations, which are simple to solve, but which 
require an uneconomically large number of time steps of limited size; and 

(2) implicit difference equations, which do not limit the time step, but 

which require at each time step the solution by iteration of a large set 

of simultaneous equations. Various finite differencing methods have 

been presented and are discussed in detail by Rosenberg [15] and Nogotov 

[16]. Convergence and stability are two Important considerations In 

using a finite difference scheme. The accuracy of any finite difference 

scheme can be controlled by choosing the space and/or time interval as 

* 

small as possible at the expense of Increased computation time in 
» 

solving the resulting increased number of algebraic equations. A finite 



difference representation of the de-icer pad Is shown in figure 2. 

Stallabrass [5] and Gent and Cansdale [17] employed the explicit 

forward time, central space finite difference technique in their de-icer 

investigation. In this numerical method, the temperature at each node 

is calculated directly from the value of the temperature at the previous 

time step. This method is stable only if certain criteria relating the 

time step and the grid size are satisfied. The criteria requires 

that the ratio of aAt to (Ax) 2 must remain less than or equal to 1/2, 

2 

i.e., the size of csAt must be of the same order of magnitude as (Ax) 
for the solution to be stable. A time step of 0.001 sec. or smaller is 
normally used so that the stability criterion at the most sensitive of 
the nodal points is satisfied. The excessive number of calculations 
needed because of this small time step can, however, cause an accumulation 
of truncation and round-off error. The truncation error for the forward 
differencing scheme is first order for the time differential and second 
order for the space differential. Stallabrass [5] used the above method 
to solve a square grid (i.e.. Ax “ Ay) two-dimensional model with a 
point source or a finite thickness heater with constant heat output. 

Baliga [4], Marano [6] and the present study have employed the 
Crank-Nicolson implicit finite difference scheme in the electrothermal 
de-icer pad analysis. This numerical scheme is unconditionally stable, 
and consequently permits an independent choice of the time and space 
parameters of the grid. The truncation error of this method is second 
order for both the tlx j and space differential. This allows a time 
step of 0.1 sec. to be used. In addition, Rosenberg [15] has pointed 
out that the Crank-Nicolson implicit finite difference scheme is the 
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preferred method for obtaining numerical solutions to pig diabolic partial 
differential equations. Ballga [4] used the Gaussian Elimination method 
to solve the resulting set of linear algebraic equations for a one- 
dimensional model of the de-icer problem. In order to obtain the 
temperature profile at any time step with the effect of the phase change, 
the Gauss-Seidel method, which is also known as the Liebmann method, is 
recommended. This method is best suited to the solution of systems whose 
matrix is sparse but not tridiagonal, i.e., for systems of two- and 
three-dimensional problems. 

Phase change or moving-boundary problems are important in many 
engineering applications, and have been studied by several investigators 
over the last two decades. Most studies have been restricted to one- 
dimensional cases. The few investigations that do deal with two- 
dimensional problems are restricted to specific problems and to a limited 
range of parameters. The solution of Stefan problems is inherently 
difficult because the interface between the solid and liquid phase is 
moving as the latent heat is absorbed or released at the interface; as a 
result the location of the solid-liquid interface is not known apriori 
and must be determined as a part of the solution. Many solution approaches 
such as Integral methods [18], Successive approximation [19], Series 
solution [20J, Green's function [2]] and Perturbation techniques [22] 
have been proposed for solving one-dimensional problems with a number 
of simplifying assumptions. The exact solution of a more general phase 
change problem was discussed by F. Neumann in his lectures in the 1860's, 
but his lecture notes containing these solutions were not published 
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until 1912. Since then, many phase change problems have appeared in the 
literature, but the exact solutions are limited to a few Idealized cases 
involving semi-inf inlte or infinite regions, subject to simple boundary 
and initial conditions [7]. Because of the non-linear nature of such 
problems, the superposition principle is not applicable and each case 
must be treated separately. Since exact solutions are not available for 
these problems, approximate, semi-analytical, and numerical methods are 
used to solve phase change problems. An extensive review of these 
techniques has been presented in References [23, 24, 23, 26]. 

Numerically, the major difficulties are caused by the need to track 
a continuously moving phase change boundary through a region approximated 
by a finite number of points. Various methods have been proposed in 
Reference [23] to overcome this problem, including moving grid points, 
isotherm migration, reformulation of the problem into Integral form, etc. 

A time dependent grid spacing scheme, also requiring nodal iteration, has 
been presented by Heitz and Westwater [27] and by Murray and Landis [28]. 
Crank [29] and Ehrlich [30] employed higher order space differences for 
the phase change equation. The drawback with this method is the require- 
ment that the media be homogeneous, thus necessitating intricate programming 
in order to simulate the movement of the phase interface across element 
boundaries . 

Most of these methods use predictor-corrector techniques to locate 
the position of the solid-liquid interface. In general, such procedures 
require a large number of trial and error computations. In order to 
circumvent these difficulties, the Enthalpy method, whieh is also called 
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the method of weak solution, has been employed frequently in the recent 
literature, e.g., see [23» 31, 32, 33]. In this technique, enthalpy is 
used as a dependent variable along with the temperature, the interface is 
eliminated from consideration in the calculations, and the problem is 
made equivalent to one of nonlinear heat conduction without change of 
phase. Hence, the advantages of this method are: (1) that the problem 

to be solved is formulated in a fixed region; (2) that no modification 
of the numerical scheme is necessary in order to satisfy the requirement 
of accurately tracking the moving phase change interface; (3) that the 
temporal derivative of the temperature no longer appears in the governing 
equation; and (4) that it is not necessary to consider the solid and 
liquid regions separately. In addition, Fox [34] asserts that weak solu- 
tion techniques are by far the best of the finite difference techniques 
for any significant problem, and especially those in more than one space 
dimension. It is shown in Reference [33] that the mathematical repre- 
sentation of the Enthalpy method is equivalent to the conventional 
formulation of the conservation equations for the solid and liquid regions 
and the movement of the solid-liquid interface. 

Kamenomostskaja [35] and Atthey [32] have applied the Enthalpy method 
to one-dimensional phase change problems using a finite difference scheme. 
The latter applied the method to a one-dimensional spot-welding problem. 
Crowley [36] applied the method to Stefan problems in more than one 
dimension. Marano [6] used the Enthalpy method in a one-dimensional de- 
icer pad analysis. 
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Bonacina et al. [37] and Goodrich [38] have proposed a heat capacity 
formulation of the Stefan problem by associating the latent heat effect 
with a finite temperature Interval about the phase change isotherm. 

Baliga [4] employed this method In a one-dimensional electrothermal de- 
icer pad analysis. Stallabrass [5] implemented the phase change process 
by holding the desired node at the melting point until enough heat had 
been transferred into the nodal volume to completely melt all of the ice 
in the volume. Both Baliga's [4] and Stallabrass' [5J techniques are 
equivalent but less rigorous than the Enthalpy method. Meyer [39j and 
Shamsunder and Sparrow [33] have used an Enthalpy technique to solve two- 
dimensional phase change problems. 

In the present study, a tra./.vient , two-dimensional solution is obtained 
for the de-icer pad by the use of finite-difference techniques which 
employ the Enthalpy method to simulate the phase change occurring in the 
ice layer. 
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III. FORMULATION OF NUMERICAL METHOD 


A. GOVERNING EQUATIONS AND BOUNDARY CONDITIONS 


The following assumptions were made in the development of a two- 
dimensional, transient, mathematical model for heat conduction In a 
composite blade with an ice layer: 

(1) The thermal physical properties of the material composing 
each layer of the blade are different, but do not depend on 
temperature; 

(2) The ambient temperature, blade air temperature and all heat 
transfer coefficients are constants; 

(3) There is perfect thermal contact between each layer; and 

(A) The density change due to melting 1 b negligible, i.e., the 

effect of the volume contraction of the ice as it melts is 
neglected, 

With the above assumptions, the mathematical formulation for the 
problem of unsteady heat conduction in a chordwise two-dimensional section 
with electrothermal heating can be represented as 


P J c pj at 


K j ( 2 
3 3x 


L > + q. 


where j stands for the material in question, as given by 
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j ■ Is blade or substrate q^ ■ 0 
j * 2s lower or inner insulation n 0 
J “ 3: heater q 3 - q 3 <t) 

j ■ 4; upper or outer insulation q^ * 0 
j » 5: abrasion shield q^ = 0 


and where 



density of the j layer; 

specific heat capacity of the J layer; 

temperature in j 1 ' 1 layer; 

thermal conductivity of the j n layer; 

rate of heat generation per unit volume in 
,th , 

j layer; 


t = time variable; and 
x,y = space coordinates. 


( 2 ) 


For the ice layer (j = 6), the governing equation for the Enthalpy 
method is: 


3H, a 3T fi a 

at - “ ax (k 6 9T" ) + ay <k 6 ajr 5 


(3) 


where 

H, = enthalpy per unit volume within the ice-water layer; 
b 

Tg “ temperature within the ice-water layer; and 
Kg ° thermal conductivity within the ice-water layer. 

The enthalpy within both phases (ice and water) can be found from 


15 


1 


equation (3). The temperature T fi can then be determined from the 
following relationship between Hg and T^: 


a C T 
d pa 6 


, T, < T 
' 6 m 


H, 


(4) 


C 1 C pl <T 6 ‘ V + "l <C P s T m + L f>' T 6 ” T m 


where 




physical properties of the solid phase; 
physical properties of the liquid phase; 
melting temperature; and 
latent heat of fusion per unit mass. 


From equation (4), 


T 6 - 


V. c - 


ps 


m 


, H, < H 
6 — sm 


,H < H, < H. 
sm 6 lm 


(H,-H. )/p,C . + T , H, > H. 
'6 lm 1 pi m ’ 6 — lm 


(5) 


with 


H 


sm 


H 


lm 


p C T 
8 ps m 

p. (C T + l f ) 
1 ps ra f' 


where H and are the enthalpies of the ice and water at the melting 
point, respectively. 
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The boundary and initial conditions are as follows: 

(i) At interior interfaces, the temperatures and heat fluxes are 
continuous, l.e., 




j ■ 1, • ■ • , S 


3T 


-K 


j 3 y, 


3T 


-K 


■j+1 3y 


j+i 


J “ 1» 5 


( 6 ) 


where I denotes the interface* 

(il) At the interior and exterior surfaces of the composite body, 

Newton's law of cooling can be used to represent the convective 
heat exchange. 

For the lower or inner boundary: 


3T 


K 


j ay. 


h bl (T j 


- T w ) . j - 1 


(7a) 


Where 1 denotes the lower or inner ambient boundary, h fel is the 
convective heat transfer coefficient at the boundary and T bl is the lower 
ambient (blade air) temperature. 

For the upper or outer boundary: 


3T 


-K 


j »y 4 


h b2 <T j 


2 “ T b2^' 


(7b) 


where 2 denotes the upper or outer ambient boundary, h b2 is the convective 
heat transfer coefficient at the boundary and is the upper ambient 
temperature. 
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(Ill) At the centerline of the heater and the gap, insulated boundary 
conditions can be used to represent the symmetry, i.e.. 


9T 

dx 

3T 

3x 


b.J - 0 

h,j - 0 


, J ■ 1, . • • , 6 
, j ■ 1, • . . , 6 


(8a) 

(8b) 


where g denotes the position of the centerline of the gap, and h 
denotes the position of the centerline of the heater. 

In addition to the above, the initial temperature distribution in 
the composite blade can be constant or a function of position, and, if 
desired, constant temperature boundary conditions can be used for the 
lower and upper surfaces of the composite blade. 

To formulate the above equations in terms of dimensionless tempera- 
ture and distance, the following definitions are made: 


8 







where 

T ref ® the reference temperature (taken to be 32 °F in this 
study); 

L = the total length of the composite slab in the y 
direction; and 

0 ^ » the thermal diffusivity of the j layer. 

Using these definitions, equations (1) through (8) take the 
following form: 

For each layer in the composite blade: 


18 


ORIGINAL PA.CX’ GS 
OF POOR QUALITY 


2 

9 6, 


30 o 3 2 0. 

3t ,2 1 -2 ,r 2 

L 9x 3y 


P j C pJ T ref 


(9) 


For the ice layer: 


3H, 

e 

3t 


r - a 30 fi a 30, 

t~ <k 6 ~ < k 6 "“)! 

L 3x 6 3x 3y 6 By 


( 10 ) 


H, 


pC 0, T c 
r a ps 6 ref 


, 6, < T 
6 ra 


p. C , T , (0, - 0 m ) + p. (C 0 m T + L.) , 0 > T 

1 pi ref 6 m 1 ps m ref f 6 m 


( 11 ) 


>H,/p C T - 
6 ps ref 


• H 6- H 


sm 


( 12 ) 


0 £ “ ( e m 
6 to 


, H < H, < H- 
sm 6 lm 


V't V T ref - f- <C P8 e B "lm 

Pi 

At the interfacial points: 


3 J+1 


30 


-K 


d 


30 j 


- -K -J±i 
*J+1 


3y 


d+i 


* J “ 1» • • • » 5 


i J “ 1> •••,5 


(13a) 

(13b) 


At the lower and upper boundaries: 
30, 


K 


j 3y 


30 


-K 


'1 - 
J ay. 


h bl 1 (6 J 


h b2 L ^ 8 j 


1 ’ 6 bl ) * J ° 1 


2 ~ 0 b2^ * J “ 6 


(14a) 


(14b) 


19 


Finally the insulated boundary conditions at the centerlines of the 
gap and the heater are; 

j - 1 6 (15a) 

8»j • 0 

J - 1, .... 6 (15b) 

h,J - 0 

Next the above equations will be expressed in finite difference 
form and arranged for numerical solution. 

B. numerical technique and crank -nicolson finite 

DIFFERENCE FORMULATION 



Formulation of a difference problem begins with the choice of a 
difference grid (a discrete set of points) , on which a discretized or 
finite difference version of the original differential equation is to 
be solved. A significant feature of this technique is associated with 
the optional choice of grid parameters for characterizing the mesh 
spacing. For this study, there are three independent variables: two 

spacial coordinates x and y, and time t. The respective grid spacings 
and time increment are Ax, Ay and At. Let E^ k stand for the value of 
a dependent variable at the point (1, k), i.e., either T, 6 or H. The 
central difference approximations of 3 E/3x and 3E/3x are: 


2 

3 E 


i.k 


3 x 2 


E i+ i j- “ 2E. k + E , . 2 

i-iiJi + o (Ax) 

(Ax) 2 


j 


( 16 ) 


IT 


'is 


3E 


i.k 


3 X 


s 


■ i ^ - k = + 0 (Ax) 2 

2 Ax 
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(17) 


These are second-order finite difference analogs. The second order 
finite difference analog for the tisae derivative, 3E/3t, is: 

A 


3E 


3t 




° k 2 

j " At -' ^ + 0 (At> 


(18) 


where the superscript A denotes the value at the previous time level and 
the superscript o denotes the value at a half-time step between the 
previous and present time levels. Equation (18) represents the Crank- 
Nicolson formulation of the time derivative. The real key to the Crank- 
Nicolson formulation is the manner of approximating 3 E/3x and 3E/9x 
without requiring the evaluation of the dependent variable at a time level 
half-way between the known and the unknown time levels. These derivatives 
are approximated by the arithmetic average of the finite difference 
analogs at the previous time level and the present time level. Equations 
(16) and (17) are approximated as: 




3x 2 


JA 


- 1/2 ( 


3 2 E 


I.k 


3x 2 


3 2 E 


3x 2 


i.k 


) 


E . .-2E, ,. + E. , , + eA . - 2E* + E 


lifcIA 


i.k 1-1. k "l+l.k 


2 (Ax) 2 


ijk 


+ 0 (Ax) 2 
(19) 


3E 


2A 


3x 


° 3E i k 

» 1/2 (— 

J 


3x 

- E, 


3E 

+ 

j 3x 


itiL 


j 


+ E* - E* 


1+1, k -i-l,k “J&jk i^ A + 0 (a - } 2 
4 Ax 


( 20 ) 
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The Crank-Nicolson finite difference equations for the present 
problem are obtained by substituting the analogs of the above derivatives 
into the governing equations and boundary conditions. 


1. Finite Difference Equations for the Composite Body 


Substituting equations (19) and (20) into equation (9) yields: 

e i,k,r e ifh..i “i+i.k.r 2e i,k.j + e i-i,k.j + e i«,k.r 28 A.j 

— ® SS? 

+ 6 l-l.k,J 8 i,k-H,j" 2 °i,k,J + 6 l,k-lJ + 8 i,k+l.j" 2 8 i»k,J + 8 i,k-l,J 


9i 

+ — s-* 1 - 


p.C ,T - 
j PJ ref 


<^> 2 


J“l» . . . , 5 


*] 


(21a) 


Solving for the unknown temperature at node (i, k) in the J th 
layer gives: 


A j { 2 tuS.\ 2 l0i +VM + Vi,k,j + Vl,k,j + e i-l,k,j 1 




2(LAx) 


te 


+ 9 , 


,+ e. A . . 4 + e, A 


[e 


,+ e 


,+ eA .. ,+ e, A 


J + B, e J 


2 (LAy ip k+lp J i|k-X f j l»k+lpj i,k-ltj J i t k f j 


q/VUy 

+ } i 2» . . . , n-1, j = 1 5 

P J C pJ T ref 


(21b) 


where 




li 
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At 


_1 

At 


<la5) 2 


1 


(LAyj ) 2 


, . V . Jj 

(LAx) 2 (LA^) 2 


and 


* the heat source per unit area, and equals q^LAy^. 

Equation (21b) Is valid for all grid points within each of the layers 
of the composite aircraft blade, except at the centerlines of the gap and 
the heater, i.e., at 1 ® 1 and 1 ® n. Finite differencing equation (15) 
gives: 


0 


0,k,j = 6 2,k.j i 9 1 6 


(22a) 


0 , . , . ° 0 , . . j « 1, . . . , 6 

n+l,k,j n~l,k,j 


(22b) 


The grid points within each of the layers of the composite body at 
the centerlines of the gap and the heater are shown in Figures 3a and 3b, 
respectively* The equations for these points can be determined by 
substituting equation (22) into equation (21b), producing the following: 
At the centerline of the gap: 


i » .1, 0 


l.k.j “J 


A, f 


J (2 0,^+2 0* J + L-t? t0 


2 (LAx) 


-,2 2,k, j 2,k,j 


2 (LAy^ ) 


- \2 l,k+l,j 


q!?/LAy, 


+ 0 +0^ +0^ 1 + B 0 ^ + ■■■J. . } 

+ e l,k-l,j + e i,k+l,j + 6 l,k-l,j + B j °l,k,j „ c T 

P j C pjref 

j = 1 5 (23a) 
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At the centerline of the heater: 

a 
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i « n, 0 , * A, { -J-s- [2 6 . , , + 20 A . , .] + =• [6 . 

n.h.J J 2 (LAk) 2 n-l,k»J 2(LAy,) 2 n,k+1 * 

J 

+ 6 , t + 0 A , . + 0 A ,,.]+B. 0 A , .} j “1, i . . t 5 
n, k-1, j n f k+l f j n|k*l|J J n*k t j 

(23b) 


2. Finite Difference Equations at Interfacial Points of the 
Composite Body 


Let (i, k) denote one of the Interfacial nodal points between the 
slab layers j and j+1, as shown in Figure 3c. Substitution of equation 
(20) into equation (13b), with the aid of analog (13a), yields: 


6 i,k,j " 6 i,k,j+l * 9 i,k,j " 0 i,k,j4 A . » j “ l ' **•» 


(24a) 


and 


0 ^ Q lQ 0 M Q U 

_ K ( i.k+l.j i.k-1.,1 i.k+l.j i.k-lJ ^ 
J /. A™ 


4 4,J 

„ . Vk+l.J+l ~ Vk-l.J+l * Vkkl.j+l 

'Vi 


- 6 


4 Ay 


i.k-l.J+1) 


J+1 


(24b) 


The temperatures 8 1>fcfl(j . 8^^. '’t.k-l.j+l a ” d 0 i?k-l,j+l 
are fictitious and can be eliminated by equation (21b) when applied to 
the layers j and j+1: 

. . - .2 

a, 


^ 2 (LAy ^ “ j 

6 i,k+l,j + 6 i,k+l,j A 4 a, { 9 i,k,j " A j [ onAZ ^2 ' u i+l,k,j 


(6, 


j J 


2 (LAx) 
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3 i-l,kJ + Vl,k,j + 9 i-l t k,j > + 2^)2 (e i,k-l,j + 6 i,k-l.j ) 


+ B 


6 q^/LAy 

i e + -■* — J t \ 

J i.k.j + 

P J°pJ T ref 


(25) 


a 2(LAy j-n ) { a i+i 

i,k-l,j+l i,k-l,J+l j+l ^ 2 (LAx) 2 


(0 t+l,k,j+l + 6 i-l,k,j+l + 6 i+l f k,j+l * 9 i-l,k,j+l ) 




+ Vl* 




A + J U+V" 7 J+1, . 


(26) 


°j+l C pj+l T ref 


Eliminating the fictitious temperatures by combining equations 
(24a), (25) and (26) yields: 


0 


i,k.J 


,2 


n (B? i >2 , Vi 6 y ~i ( *Vi )2 . 
2 1 *■ - 2 r - -2^ 


2(L& yi ) fc K 1+1 Ay^ 2(LAy 1+1 )" (Ax)" Ay 1+1 (Ax) 

————— 4 . f* . “ 


A j “j Kj Ay j+1 A j+ 1 “ j+ l 


(e i+l,k,j + 9 i-l,k,j + 9 i+l,k,j * 9 i-l,k,j I + 2 [0 i,k-l,j + 0 i,k-l,j 


2 K i+1 Ay i 

+ — iii — i. [e 
K j Ay j+i 


i,k+l,j+l + 6 i,k+l,j+l 1 + ' [ a 


2(LAy ) 2 

.. T . 


2(LAy )" K Ay 

+ i±K J - t 1 — l B ] e A + 

a K Ay j+ X ifkfi a 

j+l j y j+l 


j 


2(1.4^) 2 qV/LAy^ 


P,C .T f 
j pj ref 
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a j+l p j+l C pj+1 T ref 


i°2. 


i • i I 


n-1, J-l, ..., 


5 

(27) 


Equation (27) is applicable for the interfacial points within the 
composite blade, except at the centerlines of the gap and the heater. 
Applicable equations for the interfacial points at the centerlines of the 
gap and the heater, as shown in Figure 3c with i => 1, and i = n, can be 
determined by substituting equation (22) into equation (27): 

At the centerline of the gap: 




2tLA?j) 2 | K m 2(LAy M y 


{ 2 ( 


(Ay.) 2 

1 &- 


a j “j 


K j ay j+i A i+i “j+i 


+ lB 

(Ax) Kj Ay j+1 


2 ,k, j + 9 2,k,j J + 2 l ®l f k-l,J + Vk-i.j 1 


K 1+l Ay 1 A 2 (LAy .) 2 

+ 2 - [0 l,k+l,j+l * 0 l,k+l,j+l 1 + [ a B j 


K j Ay j+1 

2'LAy ) 2 K Ay 
+ ^ i- B . , . J 0 u 


°j+i K j Ay j+1 


j+l J 'l,k,j 
“j+l K j Ay j+1 P j+1 C pj+1 T ref 


j 

2 0-A qV/LAy^ 


j P j C P j T ref 


(27a) 


At the centerline of the heater: 


i=n. 


e 


n,k,j 


2(L A y,) Z 

— - — * — + 
Aj 0j 


«J±1 Ay 1 

“j+l A j+1 K j ^j+l 


( 2 


(Ay.) 2 

[ 

(Ax) Z 
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(^4 + l) ^<+1 A A 

+ — m iti — 1] [e + e A . .] + 2 [e . . , + 0 , . , J 

n*l,k ( j n^k— 1»J n,k—l»j 


(4x) K j 4y j + i 

2K^,iy 


2(LAy.) 


+ ■ 3 (0 , . . + 0 \ , J + [- — B 4 + 

k a Z n-l,k,J n-l,k,j ot. j a 

J fly j+1 J 


2 (LA y j+1 ) K j4.! ^ 

j+i Kj Ay j+1 


4 W?/ 2(L4? H1 f K h1 ^ 

* I- I ' ^ ^ V 


B j+1^ 0 h,k,j a 


i p j C pj T ref “j+1 K j Ay j+1 p j+l C pj+l T ref 


(27b) 


3. Finite Difference Equations for the Substrate-Ambient 
Interface 


The grid under consideration is shown in Figure 4a. At this 
surface j = 1, and (i t k) is taken to be the interfacial node. The 
finite difference representation of equation (14a) is: 

0 -0 +0 A -0 A 0 +0 A 

K (— izlii •*■■. .-jjjjJL .-- ijOg K u ^ (..jLikilLm ijJLti. . q ) 


4 Ay 


(28) 


Temperatures 0 ^ and 0 i o 1 are f ictit * ous an d must be eliminated. 
Applying equation (21b) with k ® 1, j = 1 and q^ « 0, gives: 


9 i,l,l “ A 1 { 




+ 0 , 


+ 0, 


+ 0, 


2 (LAx) 2 1+1 » 1 ' 1 i-l.-l.l i+1 » 1 * 1 ”1-1.1, 1 J 


i 0 4 , 1 +e.m + 0/0 , + 9 A 0 A 1 } (29) 


2(LAy x ) 2 i ’ 2 * 1 i* 0 * 1 i> 2 ,l i»0,l J 1 i. 1 . 1 


The fictitious temperatures are eliminated by substituting equation 
(29) into equation (28): 
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1,1,1 2(LAy 1 ) 2 Zh^LA^ (Ax) 


A 1 °l 


K, 


Uy,) 2 

( — kr- (e 


+ 0 


+ e, 


+ 9, 


-.2 ' v i+l,l,l i-1,1,1 i+1,1,1 i-1,1,1 




+ 2(e i. 2f l + 6 i,2,l> + 1 


2(LAy 1 > 2 


B 1‘ 


2 h ul LAy 


bl**"*!, q a 


K, 


1 e M,i + 


4 h bl L4y l 

K, 


e fal } i « 2 - 1 


(30) 


Again, the Interfacial points for the substrate-ambient boundary at 
the centerlines of the gap and the heater are shown in Figure 4a with 
i a l and i = n» respectively. The equations for these points can be 
determined by substituting equation (22) into equation (30): 

At the centerline of the gap: 


i°l, 6 


(Ay x ) 2 * 

{2 —At [e, , , + e, , J + 


1,1,1 2(LAy 1 ) 2 2h bl LAy x (Ax) 2 2,1,1 2,1,1 

A i “i + ~ 


x2 


A -'"“J'l' ““bl 

2 ">1,2,1 + 6 l l2 ,l J + »’ * 


2 ‘ LAy l> . 2h bl LAy l, „ A . 4h bl LAy l 




1 


Kl 

(31a) 


At the centerline of the heater: 

1 


(A^) 2 


i®n, 6 , . = ? (2 — ■ “ r * 9 - ■ [6 - , , + 0 . . i ] 

n * X 1 1 * /— . — \ 2 « * • i" //. 2 n *1 • 1 ^ 1 n™l> ^ 

,l 2 (LA y^) 2h bl LAy l 


a i °i 


K, 


2 (LAy,) 2 2hLAy 

+ 2 ">»,2.1 + + B 1 - “I q— 1 "n.1.1 


4h bi 1 - A1, i e , 
K, bl > 


(31b) 
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4, Finite Difference Equations for the Abrasion Shield- 
Ambient Interface 

When the ice layer is not present or when the ice is shed, the grid 
under consideration is shown in Figure 4b. At this surface, j = 5, and 
(i, m) is taken to be the interfacial node. Finite differencing equa- 
tion (14b) yields: 

„ . e i.m»l. 5 ' 9 l.m-1.5 + 6 l,ni+ 1.5 ' 

5 ( 4 Ay 5 

. h b2 l - v> <“> 


Tempetatures 6^, and 0^, are fictitious and ara eliminated 
by the same procedure as used for the previous substrate-ambient inter- 
face boundary: 

1 


i,m,5 2 (LAy 5 ) 2 2h b2 LAy 5 (Ax) 


(Ay 5 )^ 

T~ -,2 l 9 i+l,m,5 + 0 i-l,m,5 


A 5 a 5 




+ e iil,m,5 + °i-l,n,,5 1 + 2 16 l,m-1.5 + “i.fl.S 1 

«h b 2 L dy 5 


2(LAy ) 2 

+ 1 Zr~ ®5 ’ 


2 h b2 L j y 5 ] 0 A + 

K c J B i,m,5 * K k 


6 b 2 } 


i=2, 


. n-1 
(33) 


At the centerline of the gap: 

1 f . (Ay 5 ) 2 

iol 0 « — 5 — {2 _ 2 

* l,m,5 2(LAP 5 ) Zh^LAyg (Ax) 

-r^- + — kT - 


[0 + 
1 2,m,5 


0 2,m,5 1 
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2 (LAy.) 2 2h. 2 LAy. » 

+ 2l6 l,m-l,5 + ^.m-l.S 1 + 1 a. B 5 Kl 1 0 l,m,5 


, 4h b2 L&y 5 Q . 
+ K 5 0 b2 > 


(34a) 


At the centerline of the heater: 

i < s y 5 ) 2 4 

i*n, 6 _ q ° 5 { 2 r- [0 . , + 6 t m J 

n,m ’ 5 2(Ay 5 ) 2 2h b2 LAy 3 (Ax) 2 * n-l,ra f 5 


Ag Og Kj 


v2 


2(LAy.) 2h LAy 

+ 2 [0 , , + e . ,] + t *— b. S£ — £j e . 

n,m-l, 5 n,m-l,5 °5 5 K 5 n,m,5 


. 4h b2 u y 5 . , 

+ » / 


K 5 b2 


(34b) 


3. Finite Difference Equations for the Ice Layer 

As mentioned before, two equations are needed in the formulation of 

the Enthalpy method within the ice-water layer (j = 6) , one to calculate 
enthilpy and one to calculate temperature at each nodal point in the ice 

layer. Substitution of analogs (18) and (19) into equation (10) gives: 


- ■i.k.e . K * !lf£ 6 


At 


2L 


f “i+l.k.6 ~ 2e i.k,6 + 6 j-l.k.6 

1 - 2 

<Axr 


* S i+l.k.6 “ 29 i.k.6 + 9 i-l.k.6 . S i.k+1.6 " 20 i.k.6 + Vk-1,6 + 0 i 

+ - 2 

(Ay 6 ) 2 


- 26. A . , + e. A 




(35) 


Equation (35) can be solved explicitly for H. , , by using equation 
(12) to relate 0^ ^ ^ to H i k 6* Thus, there are three sets of equations 


. k+1 „ 6 
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pertaining to the solid phase, the melting point and the liquid phase. 
Since the thermal conductivity for each node within the ice-water layer 
is a function of position, average values of the thermal conductivity 
have been used. The quantity KW^ is the average value of the thermal 
conductivity between two adjacent nodes (1-1, k) and (i, k) ; KW,, is the 
average value of the thermal conductivity between two adjacent nodes 
(i, k) and (i+1, k) ; KW^ is the average value of the thermal conductivity 
between two adjacent nodes (i, k) and (i, k+1) ; and KW^ is the average 
value of the thermal conductivity between two adjacent nodes (i, k-1) and 
(i, k). Employing this averaging approach and substituting equation (12) 
into equation (35), we obtain, upon rearrangement, the governing finite 
difference equation for each phase. 

For the solid phase ; 


i,k,6 


a + 


At 


2 " 

2L C p 

PS 8 


KW. + KW 0 KW_ + KW, . 
[— - v - 2 + 3 , ^ > -1 


<aib 2 


< 4 V 2 


r „ 4 , At T ref ( KW 2 Vl.k.6 + KW 1 9 1-X.k.6 , "i.fcfl.6 

{H i k 6 + 2 — [ i—J— t_j_ + -~i ~ — «- 

’ ’ 21 r (Ax; (Ay 6 >‘ £ 


+ kw, e, , , , kw. a a . kw A . . 

it l»k-l, 6 , 1- /a ^ _ a ^ \ _ 2 A „ A . 

(A5) 2 (Ax) 2 i+l.k,6 ; 


KW 3 A 

+ 


KW 


- e/,. t ) - 


— (e/ 


- e. 


(A - j2 v i,k+l,6 i»k,6 (A - j2 ^i,k,6 v i,k-l,6 


:)]> 


i ° 2} * i * } n - 1 


H « . . 

i,k,6 


< H 


sm 


Substitution of equation (22) into equation (36) yields: 


(36) 
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V* 


*1 


At the centerline of the gap: 


i “ 1 * H l,k,6 " U + n ,2 


At 


2LC p 
pe a 


KW. + KW_ KW» + KW, . 
-i _ —1 : + _3 S ]} -1 


(Ax) 2 


(Ay fi ) 2 


At T ^ 2KW- 0, 


KW- 0, 


{Hl \ , + — .sm ( r7 . . xm + ziiiiM i AiJ ii ° i,k 

1 ’ R ' 6 2L 2 (Ax) 2 (Ay,) 2 




(Ax) 2 


(0 


(*yj 


- e, J 


- .2 ' l,k+l,6 -l,k,6 


KW, 


_4 f a A A , 

Z.\ 2 ( l ' k ’ 6 m l.k-l,6 ) i } H l,M < H sm 


(Ay 6 ) 


At the centerline of the heater: 


A# . KW. + KW. KW_ + KW, . 

i=n - H .H- U * I - ■ _ / ^ - 3 - . — 2- i )‘ l » 

» V. <fi *> <^ 6 > 


a At T f 2KW, e „ i i, R W, 0„ , + KW, 0 

{ ^ _ n i 3 i nik*fl t 6 A n 

n,k ’ 6 2L 2 (Ax) 2 (Ay g ) 2 


2KW 1 A (0_\ . . - e A J KW, h 


, 1 n-l.k.6 n. k.6 

(Ax) 2 


3_ A 


(Ay 6 ) 


( 0 , 


“ \2 % n,k+l,6 n,k,6 


- 9. , J 


KW. A 

4 , Q A A ... 

_ 2 ( 0* It ft ” !f_1 ft) J ) H , , < H 

„ n»K,o n,K-l»o n,k,6 — am 


(Ay 6 ) 


with 0 i,k,6 " H l f k,6 /C ps p s T ref 1 ° 1 

For the melting point: 


A At T , KW. + KW KW- + KW. 

H. . . , - H. . . + P-{ - [-i J- + ~J ±_ j ft 


l»k»6 l*k,6 ‘ 2L 2 


(Ax) 2 (Ay fi ) 2 ' Wm 


x 

=IA 


(37a) 


f k^" 1 ; 6 


(37b) 
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, KW 2 6 i+l.k.6 1 KW 1 9 l.l. k .6 , KW 3 0 l.k+1.6 * ^4 W.« 
(^) 2 (Ay 6 ) 2 

»1 A A KW 2 A A A 

+ ~TT (0 i-l k 6 " # u 6 1 " k 6 " k 6 3 + ~T7 

(Ax) 1 i,k, ° 1 » k »° (Ax)* i,k » 6 i+1 » k » 6 (Ay 6 ) Z 

fe A - e A i KW ^ ,a * 6 it 

[0 i,k+l,6 0 i,k.6 3 " (A - j2 6 i,k,6 “ 0 i,k-l,6 } 


1 " 2» > • • i n-1, H < H, . , < H. 

’ era l,k,6 lm 

‘Substitution of equation (22) into the above equation gives 


At the centerline of the gap: 

At T KW + KW KW + KW 

i=1 * H i k 6 " H i k 6 + — r 1 H-’- 1 - 2 2 + — , • 4 1 0 

1,K,D 1»K,0 /a.T\ 2 /a“ \ 2 m 


(%> 


, . kw 3 , ™ 2 -<° 2 : k . 6 - ° i:w 

( & *) 2 (Ay 6 ) 2 (ax) 2 

kw„ a . . kw. a . A 

+ ^ f 0 A _ a ^ 1_ j r a ^ ft ^ I 1 

(Ay fi ) 2 1 * k+1 ’ 6 1 * k * 6 (Ay 6 ) 2 1 * k » 6 


A A 


H < H. , , < H, 
sm l,k,6 lm 


(39a) 


At the centerline of the heater: 

A At T KW + KW. KW + KW, 

l=n, H . , » H 6 . . + |S£{- [ — i — — 2 + — 3 0 

n f k,6 n,k,6 2 L 2 ( A *) 2 (Ay,) 2 ” 


, 2KW 1 Vl.k.6 , KW 3 6 n.k-fl,6 + KW 4 6 n.k-1.6 , 2KH 1 , „ A 

<^ 2 (Ay 6 ) 2 (to, 2 

fl A 1 ^3 A A KW a a 

n,Ic,6 3 + .2 ^nfk+1,6 ” e nfk,6 3 " (A - j2 [0 n »k,6 " e nfk-l,6 J} 
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H < « . , < H. 

am n,k,o lm 


(39b) 


with 9 


i,k,6 * 0 m 


X| • * • | n 


For the liquid phase: 
H 


. KW + KW KW + KW, . 

1 k 6 ^ 2* ‘ • 2 ^ J y ) ) X 

l l K l° It . / «.T\ * /A* t* 


2* 1 - 2 
2h C plPl (Ax)* 


<*6> 


A At T . KW. + KW- KW, + KW. C 6 L. . . 

{ * i k 6 + — r * t <— r - 2 — + — — A (.fiftja + ' V 

2L <**> <^6> C pl C P l T ref 


+ KW 2 Vl,k.6 + KW 1 e i>l,k t 6 t KW 3 0 l.k-fl.6 ^ KW 4 6 l.k-1.6 


(Ax) 2 * * 


(Ay 6 ) 2 


KW, 


KW, 


+ 1- (Q A . 0 “ , --2 (Q A a A 

(Ax) 2 ( ^) 2 ^ l.k,6 " 0 l+l,k,6' 


KW“ 

+ (6 


(Ay 6 ) 2 i,k+l, 6 i,k,6 


KW 

_ a * > . A A A ... 

W U c) _ O N^4 f, £ ~ 1. t jc)]j 


(Ay^) 2 ^*k»6 i,k-l,6 


1 « 2, .... n - 1 H i,k,6 ~ H lm 


(40) 


By substituting equation (22) into equation (39), the following 

are obtained: 

At the centerline of the gap: 


At 1 9 

i=1 » H 1 k 6 " {1 + ~2^ 1 — ■ ~ 2 * + 


2L“C plPl (Ax)' 


KW, + KW„ KW, + KW. , 
3 4 ]} -l 


(Ay 6 ) Z 


AST KW, + KW- KW, + KW. C « 

w ru A + - re * r / I 2 , 3 4 \ , os p ti 

X{H l,k,6 + - 2 ((- _ 2 +— T 

2L (&X) <6y 6> C pl C pl T ref 


-) (-JS3 + 


V 
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2KW„ 6 


KW„ 6 


(Ax) 
A 


<^6> 


♦ zme . ™3 v i,k+i.6 4 kw a 9 i.k- M 2KH , 

/aH\2 , »2 


2KW n 


(0 


(AX) 


2 2.k,6 %k,6 


" A, „ J 


KW„ 

+ (6 * 


KW, 


A A 

(A? a ) 2 V ^' k+1 ' 6 ’ “ (A y )2 (0 l,k,6 * 9 l!k-l,6 )J} 


Vk.6 - H Xm 


(41a) 


At the centerline of the heater: 


4 „ , At *^1 + KW 9 kw , + KW y . 

i " n » H « k ft * U + ■■ 7 ■ » " [— * 2. + —3 4 n -l 

n ' k » 6 0.2* 1 „ *-*2 + - . 2 I' 


2L C pl ^ (AS) 


(Ay 6 )' 


/„ A AtT K V W 9 + KW, C 6 L* 

* ( Vm + -dr 1 “-hr 1 * - 1 . 2 *) <-g-a+ 

2L 4 (Ax) ' 2 


<4 V“ C P 1 C pl T ref 


- 6 ) 
m 


2KW, e 


* * W « ^n.k-l.fi ««,- , . . 

«*) (ay,) 2 w ;,2 1 n-i,k.6 ' Vk.e' 


+ kw g e 

£ 


2KWj 

(Ax) - 


KW, A A KW 

+ (Q A _ 6 A V __ 4 yg A A 

(Ay 6 ) 2 n»k+l,6 n,k,6 )2 (0 n,k,6 " 6 n,k-l,6 


)]) 


H n»k,6 - H lm 


(41b) 


with 


9i ’ M " I ref i pl»l ’ C <CpS % * ' f te7’ + *” 4 ‘ 1 


a. Abrasion -Shield-Ice Interface 

Let (i, k) be one of the abrasion shield- Ice Interface nodes as shown 
in Figure 3c (j , ■ 5). Substitution of analog equation (20) into equation 
(13b) with the aid of analog equation (13a) and with j « 5 yields: 
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e 


i.k,S 


0 i,k ,6 ’ 6 i»k,S " 0 


i»k ,6 


and 


, K ft.Wl.S ~ Vk-l.S * 6 l.W1.3 * “l.k-l.S , 

5 4 iy 5 

. . K A.ktl.6 I Vk-1.6 * e i\+1.6* e itk.-l,6 . 

6 4 4 ? 6 


(42a) 


(42b) 


Equation (21b) with J » S la used to eliminate the fictitious 
temperatures 6 ^ and 0^^+^ 5 * and equation (35) le used to 

eliminate to fictitious temperature 0^ g and 0^^^ g* Combining 
the resulting equation with equations (12) and (42a) yields the following 
equations for the solid phase, the melting point and the liquid phase. 


For the solid phase: 
H 


K S $5 , At + KH 2 . 2KW 3 


1,11,6 VeVs'S ’ 2 c ps p / ‘ «*>* ' < 4 ? 6> 2 


■] + l }” 1 X 


AtT 

<H i 6 k ,6 + jJ 6 * lte i+l,k,5 + S l-l,k,5 + 6 l+l,k,5 + ^-l.k.S 1 


V*5 

Ayg(Ax ) 2 
KW? 


KW 0 0,.. , , + KW. 0. . . , 2K. A KW, 

4 . 2 1 + 1 , k, 6 1 . 1 - 1 , k,. 6 . (Q +0 A > + _J, 

~ _ 2 - - '1 k-1 5 i k-1 S' 

(A x ) £ A« Av X » K f Av 


4y 5 4y 6 


(Ax)‘ 


(0 


A KW “ A A 2KW- 0 4 , £ 

_a 4 « 2 4 „ A . 3 ljk+1,6 , 

"j I, i) -9 1. £ ®4J.1 lr eJ + _ 2 + 


1 - 1 , k , 6 v i t k t 6' (i -j 2 v i,k ,6 1 + 1 , k , 6 


2 KW„ 


(04 


2L 2 Ay B K A 
A \ j. 0 ] } 


, A - .2 '°i,k+l ,6 " ®i,k, 6 ) + 

(Ay g ) a 5 Ayg 


i,k,5 J 


(4y ft )‘ 


102 “ 1 H l,k ,6 - H sm 


(43) 
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Substituting equation (22) into equation (43) gives 


At the centerline of the gap: 


lnl » H l,k,6 " { J 


At Ay 5 2KW 2 2K« 3 . 


V. 1 »*Vs T + <4 V 2 


•] + 1} X 


AtT ref 2K 5 Ay S A tIV "2 v 2 k 6 

A — p. (e + e \ .) + — 

{ H l,k,6 + 2L Z Ay 6 (Ax) 2 2,k * 5 2 » k * 5 (Ax) 2 


2KW„ 6, 


2K 


2KW„0 


2KW„ 


— <e +e A ) + ZX AaMLd + ZX. (e * _ 0 a ) 

Ay 5 Ay 6 1,k “ 1>S 1,k ~ 1#5 (Ay ft ) 2 (Ax) 2 2,k * 6 1,k ’ 6 


2KW 3 A 

+ i~- (0 a 


4 . . ^5 B 5 K 5 . 4 


- J + 


tty’,) 2 1 ' k ‘ 6 ' a 5 4 y 6 


Vm > > 


H, . , < H 
l,k, 6 — sm 


(44a) 


At the centerline of the heater: 


At K- Ay 

i“ n » H - {- r 

tlfiCfO r* . a*, a n « _ » * 


2KW 1 2KW* 1 

[— ry + ■ -r-? 3 +1 > « 

/A... 


IH “ , + * 

n,k # 6 y%» 2 


C ps p B Ay 6 A 5 a 5 2C pe P s L (AX) (A V 
AtT ref 2K 5 Ay 5 4 2KW ’ 9 


(0 + 0 A ) + 

2L“ Ayg (Ax) 2 ‘ n ’ 1 ’ k ’ 5 Vl > k ’ 5) (Ax) 2 


I* 


2K. 


2KW, 0 


2KW, 


A y 5 A y 6 


4. fl A v . 3 n.k+1,6 “™1 A 

( n.k-1,5 + 0 n,k-l,5 ) + } 2 + ^2 ( Vl,k,6 


A,„ A 


2 . - 


- 0 


. 2KW- (6„ . - 0 , ,) ZlThy- B.K, A 

A . 3 n.k+1,6 n,k,6 -^5 5 5 „ A 

1 f — - o ^ 


n,k,6 


(4? 6 > 2 


“5 iy 6 


n,k,5 


]} 


H , . < H 

n, k, 6 — sm 


(44b) 


with * e i,k,6 " H i,k,6 /C ps p s T ref 1 B 1 » 


n 


VS, 
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For the melting point: 
AtT, 
2L 
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2 .- 


H - ii A x _ tT ref , r 2L Ay 5 K 5 2KW i 2KW ? 

H i,k,6 H 1.M + -TT- c— - — +7 ^ + _a 


'3 S iy 6 (Ax > < iy 6 > 


2 J 9 m 


+ (e 


Ke Ay. KW„ 0, 


X A . „ a .A aw, o 

i+l.k,5 9 i-l,k,S + e i+l,k t 5 + •l.l,k.3> + J - #b k f* 


Ayg (Ax) (Ax) 4 


+ KW, 0 . 


L-HA 6 2K 5 a 

‘V- 1.3 + 9 M-1.3> 

5 0 


KW “ . 

+ ( 0 A 


(4j) 2 - e i.k. 


KW. 


2KW„6 , 


(•. a . . - e * ) x rn.kH.6 , 2kw 3 

' i.k.6 °<+1 .tf.fi’ + *— + 


(Ax) 2 *» k » 6 i+l,k,6 

* ) t 21 ^S °5 K s A , 

1,k ’ 6 “5 Ay 6 1,k ’ 5 } 


(Ay fi ) 




(Ay ) 2 ' i*k+l,6 

6 


i * 2| < 1 1 a n * 1 H < u < u * » 

sm i,k,6 H lm <*5> 

Substitution of equation (22) into the above equation yields 

At the centerline of the gap: 

1 ■ l ’ \M ■ C.e * ^ «- ^ * % * =4, 

2L a 5 A 5 Ay 6 (Ax) 2 (Ay g ) 2 m 


+ 2^(6 +0 A , 2KW 2 6 2 .k. 6 , 

Ay, (Ax) 2 2 » k,S 2,k,S .-z 


2K, 


<0, 


6'“"' ' ' (Ax)* Ay s Ay fi ' 


2KW, 


1 * k “ 1 » 5) (AZ> 2 {Q 2M " e i!k,6 ) + 


A . - 2KW 3 “ A 


(Ax) 


(Ay,) 2 ^ 1 » k+ l» 6 °l*k»6 


- e “ ,} 


+ X 2l2 *3 B s K 5 g A 

<4^) « 5 tS 6 1,M 


H < H. , < « 
sm l,k,6 lm 


(46a) 


At the centerline of the heater: 


±f " n * " Ck .6 + <" t 
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2L 2 Ay K, 2KW. 2KW, 

-+—h + -r 1 r 1 6 „ 


+e . “i 

iy 6 (to ) 2 n ' 1 ’ 1< ' 5 "-l.k.S (to ) 2 Ay, ^ < 9 n,k-l,! 


+ 0 


'6 

A 


2KW, 


) + — ’1 / g A A . 2 ^3 

n * k - 1 » s) ( ta)2 <0 n-l,k,6 “ 0 n,k,6 ) + 


(6 


(Ay ) 2 '‘“n,k+l,6 " 0 n,k,6* 
6 


. 2kw 3 Wi , 21 -~ 

(Ay A ) 2 a_ Av n »k»5 


,2 - 


“5 iy 6 


H sm < H n,k,6 < H lm 


with 0 


i,k, 6 m 


1 ^ ° * • * > n 


(46b) 


For the liquid phase: 


H - ( “ K S ^5 , 4t KW 1+KW? 2KW 

i,k,6 ~Z 7” + ~ l TT^ + + 1 } 1 X 

C pl p l Ay 6 a 5 A 5 2L C pl p l < Ax > (Ay fi ) 2 


A 4tT . 

(h/ - ref 


l,k, 6 


KW„ 6 , 


[(6, 


2L 2 lV ° i+1 » k »5 + 0 i-l,k,5 + 0 i+l,k,5 + 0 i-l,k,5 ) 7-7-% 

Ay fi (Ax) 

+ KW. 0. 


V J i 


“"2 “1+1 k (i ™i 5 4 -1 1, c 2K C 
+ - i l . l+l »M yl_ A-l.k,6 + __J_ (e + 0 A 

Ay^Ay^ i*k-l,5 i,k-l,5 


(Ax)*' 


KW 


1 (0 A 


KW„ 


(to) 2 '” t - 1 - k ’ 6 ' 9 7.6> - ~1 < 6 7,6 - \l w -> * -rjr 


2ICW, 


2KW- e 


2 - 


- *i k 6 > + ' . 2 

1,k ’ 6 (a5 6 ) 2 


3 0 i,k+l,6 r 2L Ay 5 K 3 K W x + KW 2 2KW, 

a * 4= + “".7r + 7r^l 


5 A 5 y 6 ( A *> 


(Ay 6 > 2 
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2 


m 


1 L f 2L Ay. B K A 

<C„ 9„ +s-£-)] + 5-5-5. 6 4 


'pi 


ps m T 


ref 


“5 % 


W> 


1-2 n - 1 H, > H, 

i.k,6 - lm 

Substituting equation (22) into equation (47) gives 

At the centerline of the gap: 

At K~ Ay. 

i a 1 || G3 f 5 5 At 

* l,k,6 ” 7 ~ + 


(47) 


2KW, 2KW. . 

— - + — 2J+ 1 r 1 


C pl p l Ay 6°5 A 5 2L C pl P 1 iL * )2 ( ^6 )2 

X Vk 6 + l*~ rf <B 2 k 5 + 9 2 A k 5> + — 02>U>6 

’ ' 2L 2 Ay (Ax) 2 2,k * 5 2,k,5 fA^i 2 


(Ax)' 
A 


+ — ^ — (e + e 4 ) + 2KW 3 e i.M-1.6 2KW 2 4 , 4 4 . 

45 5 4y 6 W* > + UJ 6 ,2 V « I**' 

, 2KW 3 A _ A A 2l2 A 5q K c 2KW 0 2KW 0 

7^ < We ' *1M> ' t— - 

(A V a 5 a 5 <a*> <Ay,r 


5 5 7 6 

1 * L f 2l2 ^ B. K, A 

( 9 m --T- «ps e m + W* r - 5 - 1 “A 5 1 ) 

F ref a- Ay. 

5 J 6 


'Pi 


B, . , > H. 

l,k,6 — lm 

At the centerline of the heater: 

At K «? Ay, 

i=n, H « (t ~ ^ + ££- 

n.k,6 x ; _ A - .. , + „.2 * 


(48a) 


2KW 2KW. 

l-ry + ™~rl + i > 

/A.. \ ^ 


-1 


C pl p l Ay 6 “5 A 5 2L C pl P 1 (AX) <Ay 6 ) 


r A AtT ref 2K 5 Ay 5 A 2KW i 6 , , , 

x {H , , + [ - i (0 + 6 A ) + 1 n-l,k.6 

n * k ’ 6 2L 2 Ayg (Ax) 2 n “ 1 » k > 5 n-l,k,5 /A - X 2 


(Ax) 
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2K, A 2KW, 6 ,,, , 2KW- 

.5 _ ,a . 0 A x , 3 n. k+1.6 , 1 

Ay 5 Ay 6 n » k - x ' 5 n,k-l,5 (a^) 2 (Ax) 2 


A A 2 ^3 A a Ay- K„ 

n-X,k,6 0 n,k,6 + ^-^2 (6 n» k +1.6 “ Vk.b* “ a^A^Ay” 


2KW 2KW- . * L. 2L 2 Ay-B.K- . 

+ (Ax) 2 + (Ay^) [9m ‘ C pl (Cps 8m + + “5 \ 6 "- k .5 1} 


H , , > H. 

n,k,6 — lm 


(48b) 


wlth 0 i?k,6 


^ l. k »6 

%l p l T ref 


r ( °PB 6 m + + 6 m' iBl “ 

C pl r6f 


b. Ice-Ambient Interface 


Let (1, m) be one of the Ice-ambient Interface nodes as shown in 
Figure 4c. Substitution of analog equation (20) Into equation (14b) 


yields: 


0 — 0 +0^ ^ qA 

-K ( i,m+l,6 i.m-1,6 i.m+1,6 " i.m-1,6 . 

Am ' 


9 i m fi + 6 7 m fi 

. t / i|Di|D i.m, b a \ 

\i l < 2 - V 


Temperatures 0^ ^ and 0 A m+1 ^ are fictitious and are eliminated 

using equation (35) and employing the average thermal conductivity approach. 
Combining the resulting equation with equation (12) gives 
For the solid phase: 


At KW + KW. 2KW, 2h L . , 

H i,m»6 « + l-^T^ + ~—T + -^"1 1 

2L C ps p s < 4 *) 2 <% 6 > ^6 
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"C* + , KW 2 We * a + 2I % 


2L 


(Ax) 2 


(Ay 6 ) 2 


KW 1 6 A KW A 

(Ax) 2 ( i ~ 1 » m » 6 " “-^2 (0 i t m,6 “ 0 i+l t m,6 ) " 


2I<W, 


Uy 6 r 


2h b2 L 


<9 1. 0 ,6 - 8 i.m-1.6> + ~ <* •„ - a ^ m . 6 >]> 


^ * 2| ■ • * | n - 1 


» H 4 m * < H 

l,m,6 am 


(50) 


Substitution of equation (22) into the above equation yields 
At the centerline of the gap: 


i ° 1, H. , « {1+ — M. 

l»m ,6 «. 2 : 


t i!X^ + 2 ^, r i 

(AJ) 2 (Ay,) 2 Ay* 


2L C pe p s (Ax > ( ^ 6 ) 


2KW, 9. 


v o _ aivw. w, , , 2 KW_ 

£ L 2,m t 6 + 4 l.m-1.6 2 

(a5)‘ 

A 


{H A + — r | , f / KWn 9 ' 

9T 2 . - 2 . . 2 ' - 2 
2L (Ax) (Ay fi ) /A ” V ' 


2KW, 


2, W: 


<6 2,m.6 - 9 l,m,« ) - ^2 <9 l,m,6 ' e i,«-l,6> + *>]> 


4 V 


H, m , < H 

l,m f 6 — sm 

At the centerline of the heater: 


g2 l,m t 6' 
(51a) 


1 ■ H m ° <1 + —s 
n,m,6 2 


At 


2KW 1 x 2KU 4 2h b2 L 

l _ O + ' _ r, + — 


A AtT , 

{H A . + Sg 

n,m,6 2j 2 


2L ‘ c pa p s (A *> 2 < A y6 )2 A y £ 

2KW 1 0 „ 1 _ A 2 KW. 6 , , 

[ 1 n-l.m.o + 4 n.m~l r 6 

( A x) 2 (Ay,) 2 


] J' 1 x 
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2 KM 


L (0 A 


2KU 


- e ‘ „ .) - 


^(0 A , -e A 


(Ax) 


jj2 n~l|Hi|6 n>ui|6 (Ay ^»®i6 n|iu*l|6 


) + 


^2- l (2 0 - e 4 >» 

A ? 6 S 2 -.-.6 


H , < H 

n,ra»6 — sm 


(51b) 


For the melting points 


H 


H, 


AtT , KM- + KW_ 2KW. 2h L 

re£ <[— — S-T+ -S-] 0. 


i,m,6 - "i.m.6 2I _2 “ (iJ) 2 <4$^ 2 4^ ' “ 

KW„ 0,., _ , + KW, 0. . _ , 2KW/. 0. , , KW 


I- D.,1 , T RW, o. , _ , 4RW, V. 1 , A 

2 l+l.m.b 1 i-l,ro,6 A i,m-l»o 1 . /o a 

(Ait) 2 ' (A? t ) 2 (AJ) 2 < i - 1 '"’ 6 


KW- 


2KW. 


‘ 9 l!m > 6 ) + ^2 <6 i!m,6 ' S i+l,m,6 > + (fi - *2 (9 l!m,6 ‘ 9 l,m-l,6 > 


2h. y L « 

— — (2 0 - 0/ ,)} 

Ay, 8 2 1 * ,n,6 


i ■ 2 n - 1 H < H , < H 

’ sm l,m»o — lm 


(52) 


Substituting equation (22) into equation (52) gives 


At the centerline of the gap; 

A 


i « 1, H 


l,m,6 ° H l»ra»6 " ^2 ll (A -J 


AtT , 2KW„ 2KW, 2h L 

{[ — -i- + — A? + — ] 6. 


- v2 ' - * m 

(^a) ^ 


fgj JldBzLi. - ( e A - e A ,) 

(Ax) 2 (Ay 6 ) 2 (Ax) 2 2 ’ m ’ 6 1 ’ m ’ 6 
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2 KW.“ a a 2 h» n L . 

+ H <°1 n 6 " 6 1 m 1 < 2 * a ~ e , m 6 )> 

(Ay V 1 * m * 6 g 2 l,ra,6 

'6 6 


H < H, , < H, 
am 1,01,6 lm 


(53a) 


At the centerline of the heater*. 

A AtT C 

i “ n, H„ c - hA . p- 


4tT ref 

2KW- 

f r - » 

2KW. 

4 

2L 2 

li „ 9 ^ 
(Ax) 2 

(a? 6 ) 2 


2KW, e . , 2KW. 0 . , 2KW. . A 

1 n-l.m.6 _ 4 n.m-1.6 _ 1 A - 6 ^ ) 

/a-\ 2 ~ .2 ” /a - v 2 Vl,E,6 n,m,6 ; 

(Ax) (Ay,) (Ax) 


+ r^_ (6 ‘ ,* s > « e„ .» 

(Ay ) 2 n,m * 6 n,m-l,6 8 2 n,m,6 


H < H , < H, 
8m n,m,6 lm 


(53b) 


For the liquid phase: 


At KW. + KW, 2KW, 2h. , L - 

»..= i + + _ i _ + _ a _] r 1 « 

i,m,6 2L Z C Pj> (Ax) 2 (Ay fi ) Z Ay 6 


f A AtT ref ff KW l + KW 2 , 2KW 4 4h bf W0 
l,m,6 2l 2 (Ax) 2 (Ay & ) 2 Ay & 


f~ ^ (K '' 2 W . + % - 

PA 


2KW 4 8 l.m-l.t 

(4j 6 ) 2 


(A -) 2 (9 i-l,m,6 e l,m,6 J + (A - } 2 A,m,6 


A 2KW.A . 

- 8.4.1 <«. - 4 - 8, 


2h b2 L 


1+l -“- 6 + ( 1 -"> 6 ‘ 9l -"- 1 - 6) ' 7 T ~ (2 % ■ W I} 
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1-2, n-1, H l>m 6 > H lm 

Substitution of equation (22) into equation (54) yields. 


(54) 


At the centerline of the gap; 

At 


i - 1, H. - — {1 + , 

* l,ra,6 2 


2KW_ 2KW. 2h,„L . 

_4 , H , Da ^ 


l + 


{» 


a AtT c 
A ref 


((• 


2L " C pl p l ( ^6 )2 ^6 

2KW 2 2KW 4 4h fa2 L 


-] ) x 


l , ra ,6 ,.2 l ' /A -.2 ' /A ~ ,2 

2L (Ax) (Ay fi ) Ay fe 


+ -r-> t 0 m ‘ ^ 
c , 
pi 


A L f 2KW„ 0, , 2KW, 6, m . , 

, n r, i .... t \ , i 2 ,ra, 6 4 l,m-l,6 

( C ps 0 m + T te j (i5) 2 7^)* ' 


2KW, 


2 <e A 


A 2KW 4 

- aA .) .+ * 


(Ax) 2 ' 2 ’®’ 6 X ^ 6 (Ay 


(0, , - 0 “ , ,) 

- ^2 l,m,6 l,o-l, 6 


2h L A 

— 1 — (2 0 o - 6| ,)) } 

A - g 2 l,m,6 


H- „ , > H, 

l,o,6 — lm 


(55a) 


At the centerline of the heater: 

At 2KW 2KW 2h.,L , 

1 ° n » H n m 6 " {1 + ~2 ^ I f + K + i } 1 x 

’ ’ 2L : 2 C p (Ax) 2 (Ay ) 2 Ay 


(H 


AtT ref 2KW 1 2KW 4 4h. «L 

• m , l ( r + =5* + — &-) [0 - 

0,01,6 2L 2 (Ax) 2 (Ay fi ) 2 Ay m 


-i- ( c e + A-)j - — - - KW * e n, . »-i,6 
C pl PS m r6f ( ^ )2 ( ^6 )2 
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2KW.A 

_a # a U 


(AS) 


“ A a 2KW “ . 

2 n-l.n.6 * 9 o,m,«> + ^2 < e n , B , 6 - 


0 ) 


2 h b2 L 

— T" < 2 0 O - e A ,)J) 
Ay, 82 n*D,6 

0 


H , > h 

WfiD|6 lm 


The above equatione are now ready to be solved 


numerically. 


(53b) 
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IV. NUMERICAL SOLUTION BY COMPUTER IMPLEMENTATION 


A. Gauss-Seldel Point Iterative Method 


The Crank-Nlcolson finite differencing formulation used in this 
study results in a set of algebraic equations which have to be solved 
simultaneously at each time step. The matrix size of the system is 
usually so large that the inverse by any direct method requires an exces- 
sive amount of calculations, quite apart from any considerations of 
round-off error. Most important, in this study, is the fact that the 
material phase at any grid point within the ice-water layer must be 
determined during the calculations. Accordingly, the system of equations 
has to be solved by an iterative procedure. The Gauss-Seldel method was 
chosen to perform this iteration because of its desirable convergence 
properties. 

The Gauss-Seldel method is started with an initial approximation 
which is improved by successive back substitutions. In particular, a 
first approximation is used to calculate a second approximation, the 
second approximation is then used to calculate the third, and so on until 
some convergence cirterion is met. For example, equation (21b) can be 
rewritten in the following form: 


6 


N+l 


N 


N+l 


i,k,j “ A J { 2 (La 5) 2 l0 *+l.k,j + 9 i-l,k,j +0 i+l,k,j + 0 i-l,fc,J ] 
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J i.k.j D c T 

P j C pJ T ref 


( 56 ) 


If 6^1^ j denotes the last estimate which Is known for all values 

N+l N+l N+l 

of i and k, and 0. . . denotes the new value, then 0. . . . and 641,14 

will both have already been calculated before computing 4 • All of 

the other difference equations can be rewritten in the same form as 
equation ( 56 ) . These equations are successively solved over the whole 
grid system during each time step by making a series of passes until 
the difference between two successive pass values satisfies some specified 
small value, called the convergence criterion* A convergence criterion 
of 0 . 01 % was used in this study. A complete discussion of the Gauss- 
Seidel method is given in References [ 40 J and [ 41 ]. 


B. Numerical Program Algorithm 

Figure 5 shows the flow chart of the computer program. The complete 
program Is listed in the Appendix along with a sample input data file. 

The program can be used for both conventional English units or metric 
unite by modifying the Input-output formats. The program can be reduced 
to describe the one-dimensional thermal de-icer pad by setting the gap 
width equal to zero. A metric version of the program has been compiled 
and is available upon request. 
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Although the following modifications are not presently in the 
program, they could be easily incorporated: (1) time-varying ambient 

temperatures; (2) temperature - dependent heat transfer coefficients; 
(3) initial position - dependent temperature profiles; (4) Imperfect 
Interface contact at specified nodes; and (5) temperature dependent 
properties, although this would require the use of the Kirchhoff 


temperature. 


V. DISCUSSION OF RESULTS 


The computer model formulated In this study has been designed to 
investigate the effects of a number of parameters on one and two- 
dimensional de-icer performance. De-icer performance is measured by 
the de-icing time, which is the time required to raise the ice-abrasion 
shield interface temperature to 32 °F, and, if the phase change is 
considered, to melt a thin layer of ice. The remaining ice can then be 
shed as a result of the dynamic forces acting on the outer surface of 
the aircraft blade. 

In order to test the resulting program and algorithm, Figure 6 shows 
the one-dimensional comparison of results between the present study and 
those of Campbell [11] (an approximate analytical solution), Stallabrass* 
[5J and Baliga [4]. Layer thickness and layer materials correspond to a 
standard design de-icer pad and are shown in the inset of Figure 6. The 
temperature rise plotted on th6 ordinate is the difference between 32 °F 
and the initial temperature. E.g., if the initial temperature of the de- 
icer pad is -8°F, then a temperature rise of 40°F is necessary to bring 
the abrasion shield-ice interface temperature to 32 °F, which corresponds 
to a de-icing time of approximately 7 seconds. It can be seen that 
excellent agreement is achieved between the present program and the 
approximate analytical solution. The results of Stallabrass [5] and 
Baliga [4] show a slightly optimistic temperature rise beyond 5 seconds. 
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A number of parameters were Investigated using the present program* 
and* wherever possible* the results are compared with the results of 
Stallabrass [5], Ballga [4] and Marano [6] to illustrate de-icer 
performance* The one-dimensional model Is the subject of Sections A 
through F below. The phase change Is not considered In Sections A through 
H* whereas Sections 1 through L Incorporate the Enthalpy method for the 
phase change simulation. For all of the following cases* unless other- 
wise specified, values of h^ “ 1 Btu/ft 2 -hr-°F and h b2 » 10 6 Btu/ft 2 - 
hr-°F were used. Material property data for this study are presented in 
Table 1. 

A. Effect of Power Density 

Figure 7 shows the results of the effect of heater power density on 
the de-icing time for various ambient temperatures. Thu power density for 
the heater ranged from 15 to 40 Watts/in . Good agreement is obtained 
with the results of Stallabrass [5], Baliga [4] and Marano [6]. A 
slight variation between results occurs because of the more optimistic 
nature of Stallabrass' and Baliga' s results as indicated in Figure 6. 
Examination of Figure 7 reveals that the practical minimum power density 
for the de-icer is 25 Watts/in * which agrees with the observation made 
from tests both in the laboratory and on various helicopter blades [5J. 

For power densities on the order of 15 or 20 Watts/in , particularly at 
low ambient temperatures* de-icing times are quite long. 


B. Effect of Heater Element Thickness 


The heater element usually is made up either of a woven mat of wires 
and glass fibers, or of resistance ribbon. Woven mats are an order of 
magnitude thicker than resistance ribbon. Figure 8 presents a comparison 
with the results of Marano [6] for different thickness heaters. Curve 1 
represents the response of a point (zero thickness) heater, which is an 
idealization and yields the best possible results. Curve 2 represents 
the response of a heater composed of resistance ribbon, and Curve 3 
represents the response of a heater composed of a woven mat. It can be 
seen that good agreement between the two studies is obtained. In addi- 
tion, Figure 8 shows that heater thickness does not drastically affect 
de-icer performance. The maximum difference between the de-icing times 
for the various thickness heaters shown is le.ss than two seconds. 

C. Effect of Insulation Thickness Ratio and Insulation Material 

In order to obtain optimum de-icer performance, it is desired that 
a maximum amount of energy released from the heater be directed toward the 
ice layer. This can be accomplished by increasing the insulation thickness 
beneath the heater. The effect on energy transfer of increasing this 
thickness relative to the insulation thickness above the heater is shown 
in Figure 9. An insulation thickness of 0.010 inches was assumed for the 
outer epoxy/glass insulation layer for the present investigation. Figure 
9 presents a comparison of de-icing times with the results of Baliga [4] 
for the effect of varying the insulation thickness ratio: inner insulation 

thickness/outer insulation thickness. Again, good agreement between the 
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studies Is obtained. As shown in Figure 9, when the ratio is increased 
from 1 through S (Curve 1, 2, 3), the de-icing time decreases appreciably. 

It can be seen that Increasing the inner insulation from 0.020" to 0.050" 
reduces the time required for a 45°F temperature rise from nearly 8.7 
seconds to just over 6.9 seconds, a 21% improvement. It is also seen that 
increasing the insulation thickness ratio from 5 to 10 only has an 
appreciable effect for very low ambient temperatures (Curves 3 and 4). 

Figure 10 shows a comparison with Baliga's results for the effect 
on the de-icing time when the epoxy/glae” inner insulation layer having 
a thermal conductivity of 0.22 Btu/ft-hr-°F is replaced by a layer of 
Polytrifluorochloroethylene (KEL-F) of the same thickness but having a 
thermal conductivity of 0.04 Btu/ft-hr-°F. As shown in Figure 10, for a 
45°F temperature rise, the de-icing time is reduced by 40% (Curves 1 
and 3). Therefore, the use of an inner insulation material of very low ther 
mal conductivity is recommended if it is also a good electrical insulator. 

D. Effect of Substrate Material 

An interesting result was obtained when a stainless steel substrate 
was substituted for the aluminum alloy substrate. It might be expected, 
because of the lower thermal conductivity and lower thermal diffuslvity 
of the stainless steel, that the de-icing time would be reduced. However, 
the opposite was in fact found to be the case, as shown in Figure 11, 
which also contains Stallabrass' [5] results. This result is attributed 
to the higher thermal capacity per unit volume of the stainless steel. 
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E. Effect of Variable Heater Output 


A number of time-dependent heater outputs have been investigated in 
this study. In this section, a ramp function heater output is presented 
and the temperature response at all Interfacial nodes is determined. 

Figures 12a and b depict the temperature response for a point heat source 
with a ramp function output, 5 seconds on, 1 second off. These results 
are compared to a constant point heat source. In all cases the ambient 
temperature was -4°F. Both sources have an average power density of 
25 Watts/ in . As expected, the temperature at the various interfacial 
nodes drops as the heater is switched off, and rises after the heater is 
switched on. The temperature at the outer surface of the ice remains 
constant at -4°F. This is due to the large heat transfer coefficient at 
this surface. As indicated in Figures 12a and b, if heat is applied such 
that the heater is operated for a period of time until the ice-abrasion 
shield interface reaches 32^, and then is switched off until another 
layer of ice forms on the blade surface, the total energy consumption for 
de-icing can be reduced. The computer pregram developed in this study 
has been tested with constant, ramp and sinusoidal heater outputs; however, 
it can be easily modified to accommodate any other type of variable 
heater output. 

F. Effect of Initial Ice Layer Thickness 

In order to determine what effect the ice layer thickness has on 
de-icer performance, an investigation was performed in which the Initial 
ice layer thickness was varied from 0.10" to 0.50". Figure 13 shows the 
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results of this Investigation* along with the corresponding results 
obtained by Marano [6]. It can be seen that the de-icing times for the 
thicker ice layers did not differ appreciably. However, for a thin 
ice layer (0.10") with a large convection rate at the ice-ambient inter- 
face, the de-icing time increased significantly. The reason for this 
unexpected result is that any ice layer acts somewhat as an "insulator" 
and therefore slows the passage of heat, thus allowing the supplied thermal 
energy to raise the abrasion shield-ice Interface temperature. Consequently, 
for thicker initial layers of ice, this rate of temperature rise is 
relatively high, which results in shorter de-icing times. Conversely, for 
the case of very thin initial ice layers with high convection rates at. the 
outer surface, the thermal energy supplied to the abrasion shield-ice 
interface is conducted through the ice layer and is convected to the 
ambient at high enough rates that the Interface temperature rises more 
slowly, thus increasing the de-icing time. 

G. Effect of Gap Width 

Earlier attempts at analyzing transient heat flow in de-icer pads 
with only one-dimensional mathematical models proved to be inadequate 
for prediction purposes, as a 3 to 5 times greater energy requirement 
was needed in actual operation as versus that predicted. This difference 
has been attributed to the effect of two-dimensional heat flow, l.e., 
the effect of delayed heating over the gap between heater elements. The 
two-dimensional model of Figure 2 was programmed to investigate the 
effect of the width of the gap between legs of the heater elment. Gap 
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widths of 0.010”, 0.030”, 0.050” and 0.070” were Investigated. Figure 
14 employed an Insulation thickness ratio of 2, a heater of zero thick- 
ness, but having a width of 0.510”, and a 0.010” stainless steel abrasion 

shield and a 0.100" substrate. The heat transfer coefficients used were 
2 2 

hjjl a ^0 Btu/ft -hr-°F and h^ “ 200 Btu/ft -hj-°F. Figure 14 shows how 
much more time is required to achieve a given temperature rise over the 
gap than over the heater element. The results of Stallabrass [5] are also 
shown. It can be seen that to raise the temperature over the centerline 
of a 0.070” gap by 36°F requires 5.2 seconds more time than to raise the 
temperature over the centerline of the heater element by the same amount. 
The two-dimensional results of Stallabrass [5] are seen to be more 
optimistic than those obtained in the present study. Since Stallabrass 
used an explicit finite difference formulation which required small time 
and gpace Increments in order to maintain stability and achieve accuracy, 
it is possible that excessive computer time prevented sufficient testing 
to find the increment sizes below which repetitive results were obtained. 
This would account for the difference in the results of the two studies, 
particularly at large gap widths. A 0.010” gap is seen to have an almost 
negligible effect on the abrasion shield-ice Interface temperature rise. 

H. Effect of Heat Transfer Coefficient 

The magnitudes of the heat transfer coefficients, h^ and h^ of 
Figure 1, have a greater effect on de-icer performance than do any of the 
parameters previously discussed when the blade has a very thin initial 
ice layer at low ambient temperature. As shown in Table 2, to raise the 


temperature over the centerline of a 0.070" gap by 36°F requires 11.3 

2 

seconds with heat transfer coefficient values of h^ ■ 10 Btu/ft -hr-°F 

and h^ ■ 200 Btu/ft 2 -hr-°F, and 14.2 seconds with h^ ■ 1 Btu/ft Z -hr-°F 

and h „ ® 10 6 Btu/f t 2 -hr-°F. Thr initial ice layer thickness iB 0.15". 
b2 

The reason for this phenomenon is already described in the Section F, 
i.e., the heat is rapidly conducted away from the abrasion shield-ice 
interface for very thin ice layers, low ambient temperatures, and high 
heat transfer coefficients. It can be seen from Table 2 that the de-icing 
time does not change for a temperature rise of 10.8°F. However, for high 
temperature increases (lower ambient temperatures) and greater gap widths, 
the heat transfer coefficient effect is very evident. 

I. Effect of Phase Change 

In the earlier part of this study, the phase change in the ice layer 
has been neglected, i.e., no account has yet been made for the latent 
heat needed to melt the ice. The Enthalpy method is applied to simulate 
the phase change in the ice layer as described in Chapter III. Figure 
15 shows the temperature response for the substrate, heater, abrasion 
shield-ice interface and ice-ambient interface in the two-dimensional de- 
icer pad. It can be seen that the abrasion shield-ice interface over 
the center of the gap requires more time (11 sec.) to achieve a 
temperature of 32 °F than does the same interface over the center of the 
heater element (4.6 sec.). 

Figure 16 shows a more detailed comparison with the results of 
Stallabrass [5] for the effect of the phase change on the abrasion shield - 
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Ice interface temperature. This de-icer pad had a gap width of 0.050", 

2 

a power density of 25 Watts/in , an initial equilibrium temperature of 
23°F, and a stainless steel abrasion shield and substrate 0.010" thick, 
respectively. Two curves are shown in Figure 16: the Interface 

temperature above the center of the heater element and that above the 
center of the gap, both plotted as a function of time. Excellent agree- 
ment is achieved between the present study and Stallabrass’ results. 

Curves 1, 4 and 6 show the Interface temperatures would behave if the 
latent heat of the ice were not taken into account, i.e., if no phase 
change took place. These curves show that the interface temperatures 
Increase continuously. Curves 2, 3, 5 and 7 show the temperature response 
if the latent heat of the ice is taken into account. Illustrating that a 
plateau becomes evident when the abrasion shield-ice interface temperature 
reaches 32°F. This plateau remains until enough heat hue been transferred 
into the nodal volume to completely melt all of the ice in that volume, 
after which the interface temperature smoothly increases. However, this 
plateau is a consequence of the numerical solution technique, and in 
reality does not occur. As mentioned in Marano [6], the temperature 
response at the abrasion shield-ice interface and for any point in the 
ice layer has a very strong nodal dependence. If more than 90 nodes are 
used in the ice layer, the temperature response curves will not exhibit 
a plateau. However, such simulations consume a large amount of computer 
time, and the average temperatures obtained are identical with those 
resulting when a lesser number of nodes is used [6 j . It can be seen from 
Figure 16 that the phase change does not occur over the center of the 
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gap until about 4.2 seconds for the present study, as versus 0.9 seconds 
over the center of the heater. 

J. Effect of Abrasion Shield 

The abrasion shield has sometimes been called the heat diffuser for 

it will diffuse or conduct the heat laterally, thus producing greater 

uniformity of the temperature distribution at the abrasion shield-ice 

Interface. Table 3 employed the same basic de-icer pad as given in 

Figure 16, with the abrasion shield material being changed each time. 

Good agreement is shown in Table 3 as compared with the results of 

Stallabrass [5}. Case 1 is the case presented in Figure 16, i.e., a 0.010" 

stainless steel abrasion shield. The values tabulated are the de-icing 

times for the abrasion shield-ice interface to reach 32°F over the center 

of the heat' : element and over the center of the gap. Initial equilibrium 

2 

temperature in all cases was 23°F and the power density was 25 Watts/in . 
Case 2 shows that replacing the thin stainless steel abrasion shield by 
the same thickness of a more highly conductive material (nickel) has 
little effect. This is true also of a 0.010" plating of nickel on 
0.010" of stainless steel (Case 3). However, as the thickness of the 
more conductive material is increased, a great improvement is shown (compare 
Case 3 and Case 4). Case 5 shows that by locating the better thermal 
conductor beneath the poorer one, more heat flows parallel to the blade 
surface from the warmer to the colder regions. In Case 6, an additional 
insulating layer has been inserted between the two metal layers of Case 
5, Although the de-icing time has been increased, the uniformity of the 
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heat distribution has been greatly improved. However, this process of 
improving lateral heat flow should not be taken too far because the 
added material will increase the thermal capacity of the system which in 
turn will Increase the de-icing time. 

Table 3 depicts the effect of Increasing the thickness of a good 
conductor (nickel) from zero (Case 1) to 0.020" (Case 4). Further 
increase in the thickness of the nickel will yield further Improvement 
in de-icing time over the gap before passing through a minimum; then the 
de-icing time will start to increase. Clearly, narrow heater gaps 
require a thinner abrasion shield as a heat diffuser to produce greater 
uniformity of surface temperature, while wider heater gaps require a 
thicker abrasion shield. 

K. Effect of Ice Shedding 

As soon as a thin layer of ice is melted, the ice layer plus any 
water which has formed is shed from the de-icer pad due to the aero- 
dynamic and/or centrifugal forces acting on the outer surface of the 
composite blade. Once the heater has been turned off, a new ice layer 
may form. Figure 17 shows the temperature response of the de-icer pad 
for ice shedding for the first 20 second cycle. The heater was turned 
on for 10 seconds and then off for 10 seconds, and the ice was replaced 
every 20 seconds. The ice was shed when the interface over the heater 
gap reached 32°F and melted. Sharp decreases in the heater element 
temperatures are shown to occur at 6.8 seconds and at 10.0 seconds. The 
former decrease is from the ice being shed, and the latter Is from the 
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heater being turned off. The abrasion shield sur£ace temperature drops 
immediately from 32°F to -4°F when the Ice Is shed. Only the substrate 
temperature remains warm after 20 seconds. During ice shedding, the 
program time step had to be drastically reduced in order to accurately 
follow the abrupt changes in temperature which occur. 

Figure 18 shows the temperature response for the second cycle, from 
20 to 40 seconds. The temperature profiles are quite similar to those 
in Figure 17. The only difference is that, since the substrate tempera- 
ture is initially hotter because of residual heat, the de-icing time is 
decreased to 6.2 seconds. 

L. Effect of Refreezing 

From the results of Section K, it would be interesting to see the 
composite blade temperature response if the ice layer could not be shed 
and refreezing occurs. As shown in Figure 19, the heater was turned on 
for 10 seconds and then was turned off. In this case, the water at the 
abrasion shield-ice Interface begins to refreeze after 12.0 seconds, 
and has completely solidified after 19.2 seconds. 
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VI 


CONCLUSIONS AND RECOMMENDATIONS 


The two-dimensional numerical simulation model developed in this 
study for the electrothermal de-icer pad has been shown to predict accurate 
temperature profiles for several boundary conditions and thermal heat 
sources. Parametric studies using the simulation to investigate such 
practical composite blade design aspects as the effect due to different 
heating functions on the rate and depth of ice melting, the effect of 
the abrasion shield on lateral heat transfer, the effect of the phase 
change on de-icer performance, etc., have been completed. The results 
were found to agree comparatively well with previous numerical calcula- 
tions by Stallabrass [5], Ballga [4} and Marano [6]. The simulation 
contains the following improvements over Stallabrass’ [3] two-dimensional 
numerical technique for the electrothermal de-icer pad: 

(1) The Crank-Nlcolson implicit finite difference scheme was 
used Instead of the explicit forward finite difference method. 
This decreased the total number of calculations and promoted 
accuracy and stability; 

(2) The Enthalpy method was used to simulate the phase change 
in the ice layer; 

(3) The simulation can apply to shedding and refreezing of the 
ice layer. 

(4) In addition, the program can be easily modified to handle time 
and posit ion -varying ambient temperatures, temperature-dependent 
heat transfer coefficients, and imperfect contact between layers 
at specified nodes. 
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The computer program was run on the University of Toledo IBn 4341 
computer. A 20 x 29 mesh was used to represent the 6 layers of the 
composite body. Typical computing times for a total simulation time of 
20 seconds were 11 minutes without the phase change and 18 minutes 
with the. 2 >toea change. 

Baseti on the progress thus far, it is recommended that the following 
additional work be carried out: 

(1) A complete experimental study should be conducted on de-icer 
pade in order to provide an exhaustive examination of the 
assumptions uned in the simulation model and on model accuracy. 

It was originally anticipated that such experimental data 
would be available from teste conducted in the NASA Lewis 
Icing Tunnel. Unfortunately, the experimental study, origi- 
nally scheduled for Summer, 1982, has been delayed until 

summer, 1983. At that time time important analysis should 
be carried out: 

(2) In order to assess the effect of blade geometry on de-icer 
design, an analysis should be conducted using the true blade 
shape. This work is already underway at The University of Toledo, 
where a free-surface coordinate transformation technique is 
being investigated. 

(3) Available model simulations which predict the rate of ice 
deposition on blade surfaces should be coupled with the 
internal blade heat transfer programs in order to simulate 
the total process of ice formation and removal. 
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COMPLETE PROGRAM LISTING 
AND 

SAMPLE INPUT DATA PILE 
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//UOFT157S JuB (UT, 

// 06200016,19) 

// EXEC FORTXCLG,PARM a 
//FORT.SYSIN DD * 

C 
C 
C 

THIS PROGRAM CAN CALCULATE THE TEMPERATURE PROFILE IN 
A COMPOSITE SLAB WHICH HAS CONVECTIVE .CONSTANT TEMP- 
ERATURE OR MIXED BOUNDARY CONDITIONS. 


'BCD' 


NUMERICAL SIMULATION OF TWO-DIMENSIONAL HEAT TRANSFER 
IN COMPOSITE BODIES KITH PHASE CHANGE. 


THE PROGRAM CAN ALSO BE USED FOR COMPOSITE BODY PROBLEMS 
WITH CONSTANT OR VARIABLE HEAT SOURCES. 


C 

AK 

= THERMAL CONDUCTIVITIES OF LAYERS. 

C 

ALP 

= THERMAL DIFFUSIVITIES OF LAYERS. 

c 

AKL 

= THERMAL CONDUCTIVITY OF WATER . 

C 

ALPL 

= THERMAL DIFFUSIVITY OF WATER. 

C 

CA 

* CONSTANT. 

c 

CHANGE 

= FUNCTION PROGRAM TO CORRECTING THERMAL 

c 


CONDUCTIVITIES OF EACH NODE IN THE ICE LAYER 

c 


FOR PHASE CHANGE 

c 

CPS.CPL 

= SPECIFIC HEAT OF ICE AND WATER PER UNIT VOLUME. 

c 

DES .DEL 

= DENSITY OF ICE AND WATER. 

c 

DTAUI 

= INITIAL TIME STEP. 

c 

DTAUM 

« INTERMEDIATE TIME STEP. 

c 

DTAUF 

« FINAL TIME STEP. 

c 

DX 

= .'PACING BETWEEN NODES IN THE X-DIRECTION OF THE 

c 


GRID. 

c 

DY 

= SPACING BETWEEN NODES IN THE Y-DIRECTION OF THE 

c 


GRID. 

c 

EL 

= LENGTH OF EACH LAYER. 

c 

HEAD 

= HEADINGS. 

c 

HI 

= HEAT TRANSFER COEFF. AT LOWER BOUNDARY. 

c 

H2 

= HEAT TRANSFER COEFF. AT UPPER BOUNDARY . 

c 

HLAM 

= LATENT HEAT OF ICE. 

c 

IBC1 

« 1, IMPLIES TEMPERATURE IS CONSTANT AT X=l. 

c 

IBC1 

« 2, IMPLIES CONVECTIVE HEAT TRANSFER AT X=l. 

c 

IBC2 

= 1, IMPLIES TEMPERATURE IS CONSTANT AT X=2. 

c 

IBC2 

= 2, IMPLIES CONVECTIVE HEAT TRANSFER AT X=2. 

c 

ICL 

= NUMBER OF CYCLE FOR SHEDDING ICE. 

c 

ICOUNT 

= COUNTER ON TIME STEP. 

c 

IFREQ 

= NUMBER OF TIME STEPS BETWEEN SUCCESSIVE 

c 


PRINTING OF THE TEMPERATURE PROFILE. 

c 

IG 

= 1, IMPLIES PHASE CHANGE IN ICE LAYER IS NOT 

c 


CONSIDERED. 

c 

IG 

= 2, IMPLIES PHASE CHANGE IN ICE LAYER IS 

c 


. CONSIDERED, 

c 

1H 

= 1, IMPLIES NO HEAT SOURCE. 

c 

IH 

= 2, IMPLIES POINT HEAT SOURCE. 

c 

IH 

= 3, IMPLIES HEAT GENERATION WITHIN SLAB. 

c 

IHQ 

= 1, IMPLIES CONSTANT HEAT SOURCE. 

c 

IHQ 

= 2, IMPLIES STEP FUNCTION HEAT SOURCE. 
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ORIGINAL PAGE IS 
OF POOR QUALITY 


V* 


IHQ 

IHQ 

IJ 

IRF 

IRF 

ISH 

ISH 

L 

LI 

L2 

M 

MM 

N 

NISP 

NMSP 

NFSP 

NHSP 

NTSP 

NO 

NOl 

N02 

NODE 

Q 

Q2 

QHEAT2 

QHEAT3 

QHEAT4 

QV 

T 

TE 

TGI 

TG2 

TIN 

TLEN1 

TLEN2 

TLEN3 

TMP 

TOFF 

TON 

TR 

TREF 

TX1 

TX2 


a 3, IMPLIES RAMP ^UNCTION HEAT SOURCE. 

* 4, IMPLIES SINE FUNCTION HEAT SOURCE. 

» SLAB WITHIN WHICH HEAT GENERATION OCCURS, 
a 1, IMPLIES REFREEZING ICE LAYER IS NOT 

CONSIDERED. 

a 2, IMPLIES REFREEZING ICE LAYER IS CONSIDERED, 
a 1, IMPLIES SHEDDING ICE LAYER IS NOT CONSIDERED, 
a 2, IMPLIES SHEDDING ICE LAYER IS CONSIDERED, 
a NUMBER OF LAYERS IN SLAB 
a LOWER SLAB NUMBER FOR POINT HEAT SOURCE, 
a UPPER SLAB NUMBER FOR POINT HEAT SOURCE. 

= NUMBER OF NODES IN THE Y-DIRECTION OF THE GRID 
a INTERFACE NODE NUMBERS. 

= NUMBER OF NODES IN THE X-DIRECTION OF THE GRID 
a NUMBER OF TIME STEPS FOR WHICH INTIAL TIME 

STEP IS USED 

a NUMBER OF TIME STEPS FOR WHICH INTERMIDIATE 

TIME STEP IS USED. 

a NUMBER OF TIME STEPS FOR REFREEZING ICE LAYER. 

= NUMBER OF TIME STEPS FOR SHEDDING AT 2ND CYCLE, 
a NUMBER OF TIME STEPS FOR STOPPING THE PROGRAM. 

= NUMBER OF NODES IN EACH LAYER IN THE Y-DIRECTION 
OF THE GRID 

a LOWER NODE NUMBER FOR FINITE THICKNESS HEATER, 
a UPPER NODE NUMBER FOR FINITE THICKNESS HEATER, 
a NODE AT WHICH POINT HEAT SOURCE IS APPLIED, 
a POINT HEAT SOURCE WATTS/ IN*IN. 

= VOLUMETRIC HEAT SOURCE WATTS/ IN*IN*IN. 
a FUNCTION PROGRAM FOR STEP FUNCTION HEAT INPUT, 
a FUNCTION PROGRAM FOR RAMP FUNCTION HEAT INPUT. 

= FUNCTION PROGRAM FOR SINE FUNCTION HEAT INPUT. 

= INPUT FOR VARIABLE HEAT SOURCE. 

= NON-DIMENSIONAL TEMPERATURE. 

= TEMPERATURE. 

a AMBIENT TEMPERATURE AT LOWER BOUNDARY OF SLAB. 

= AMBIENT TEMPERATURE AT UPPER BOUNDARY OF SLAB. 

= INITIAL TEMPERATURE IN SLAB. 

= TOTAL LENGTH OF THE SLAB IN THE Y-DIRECTION 
OF THE GRID 
a WIDTH OF THE HEATER 
= GAP WIDTH 

= ICE MELTING TEMPERATURE, 
a OFF TIME OF STEP HEAT INPUT. 

= ON TIME OF STEP HEAT INPUT, 
a NON-DIMENSIONAL TEMPERATURE AT PREVIOUS 
TIME STEP 

a REFERENCE TEMPERATURE . 

= CONSTANT TEMPERATURE AT LOWER SLAB BOUNDARY, 
a CONSTANT TEMPERATURE AT UPPER SLAB BOUNDARY. 


DATA IN/5/, 10/6/ 

DIMENSION HEAD (42 , 80) , DIF(30 , 95 ) , THOLD ( 30., 95 ) , C0NST( 30 , 95 ) 
DIMENSION PHOLD (30 , 95 ) , TOLD (30 , 95 ) , H ( 30 , 95 ) , HO ( 30 , 95 ) 
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DIMENSION T(30,95),TR(30,95),T0(30,95),TEO(3O,95),FDIST(30,95) 
DIMENSION ALPaj.CPDaj.NOidhMMdhDYdJ.AKdhELaj.TEOO.PS) 
DOUBLE PRECISION E1,FI,G1,P1,P2 
DOUBLE T VISION A(7) ,B(7),C(7) ,D(7) ,E(7) 

COMMON /c iAl/HSMP,HLMP,AKS,AKL 

INPUT DATA 

DO 5 1=1,20 

READ(IN,7Q0)(HEAD(I,J),J=1,80) 

5 CONTINUE 

READ(IN,701)L,N,M,TLENI,TLEN2,TLEN3 

DO 10 1=1, N 
DO 10 J=1,M 
PHOLD(I, J)=0. 

TOLD(I,J)=0. 

THOLD(I , J)=0. 

CONST(I,J)=0. 

FDIST(I,J)=0. 

H(I,J)=0, 

HD(I,J)=0. 

T(I,J)=0. 

TE(I ,J)=0. 

TO(I , J)=0. 

TE0(I, J)=0. 

TR(I ,J)=0. 

DIF(I, J)=C . 

10 CONTINUE 
DO IS K=1,L 
A(K)=0. 

B(K)=0. 

C(K)=0. 

D(K)=0. 

DY(K)=0. 

E(K)=0. 

CPD(K)=0. 

ALP(K)=0. 

AK(K)=0. 

EL(K)=(j. 

DY(K)=0, 

NO(K)=0 
MM(K)=0 
IS CONTINUE 

INPUT DATA 

DO 20 1=21,24 

20 READUN,700)(HEAD(I,J),J=1,80) 

DO 25 K=1,L 

RE AD (IN, 702) NO(K) ,EL(K) ,AK(K),ALP(K) 

25 CONTINUE 

DO 30 1=25,27 

30 READ(IN,700)CHEAD(J,J),J=1,80) 
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READ(IN,703)IH,N0DE,L1,L2 
READ( IN , 704) I J ,N01 ,N02 , I HQ 
READ(IN,705)Q,TON,TOFF,QV 
DO 35 1=28,30 

35 READ(IN,700)(HEAD(I,J),J=1 ,80) 

READ(IN, 706) IBC1 , IBC2 
READ(IN,716)TX1 ,TX2 
READ(IN , 707)TG1 ,H1 ,TG2 ,H2 
READ ( IN , 708 )TIN ,TREF 
DO 40 1=31,33 

40 READ(IN, 700) (HEAD(I , J) , J=1 ,80) 

READ (IN , 709 ) IG ,HLAM , ALPL ,AKL,DEL,CPL,DES ,CPS 

READ(IN,710)TMP 

DO 45 1=34,36 

45 READ(IN,700)(HEAD(I,J),J=1,80) 

READ ( IN , 7 1 1 ) DTAUI , N I SP , DTAUM , NMSP , DTAUF 
READ(IN,712)IFREQ 

READ( IN , 713) ISH , IRF.NFSP , ICL ,NHSP ,NTSP 


PRINT THE DATA 


DO 50 1=1,20 

50 VRITE(IO, 775) (HEAD(I , J) , J=l,80) 

WRITE(I0,717)L,N,M,TLEN1,TLEN2,TLEN3 
DO 55 1=21,24 

55 WRITE(I0,775) (HEAD(I, J) ,J=1,80) 

DO 60 K=1,L 

WRITE(I0,702)N0(K) ,EL(K) ,AK(K) ,ALP(K) 

60 CONTINUE 

DO 65 1=25,27 

65 WRITE(IO,775) (HEAD(I,J) ,J=1,80) 
IF(IH.EQ.l) GO TO 70 
IFCIH.EQ.2) GO TO 75 
WRITE(I0,718)IJ,N01,N02 
GO TO 80 

70 WRITE(I0,719) 

GO TO 90 

75 WRITE(IO,720)NODE,L1,L2 

80 IF(IHQ.NE.l) GO TO 85 
WRITE CIO, 721 )Q 

85 IF ( IHQ . ME . 1 ) WRITE ( 10 , 7 22 )TON ,TOFF 
IF(IHQ.EQ,2)WRITE(I0, 724)QV 
IF(IHQ.EQ.3)WRITE(I0,726)QV 
IF(IHQ.EQ,4)WRITE(IO,727)QV 

90 CONTINUE 

DO 95 1=28,30 

95 WRITE(IO,775)(HEAD(I,J) ,J=1,80) 
IF(IBC1.EQ.2) GO TO 100 
WRITE (10, 729 )TX1 
GO TO 105 

100 WRITE(I0,730)TG1 ,H1 
105 IF(IBC2.EQ.2) GO TO 110 
WRITE (10, 732 )TX2 
GO TO 115 
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110 WRITE(I0,734)TG2,H2 
115 WRITEOO, 736)T1N,TREF 

WRITE(I0,775)(HEAD(28,J),J=1,80) 
IF(IG.EQ.1)G0 TO 120 
WRITE(I0, 738) 

WRITE(I0,775)(HEAD(31,J),J=1,80) 
WRITE(I0,775)(HEAD(30,J) ,J=1,80) 

WRITE (10, 775) (HE AD (32, J) ,J=1 ,80) 
WRITE(IO,775)(HEAD(33,J),J=l f :0) 

WRITE ( 10 , 740 )HLAM, AKL, ALPL, DEL, CPL, DES 
WRITE (10 ,750)CPS 
WRITE (10, 752)TMP 
GO TO 125 
120 WRITE(IO,7S3) 

125 WRITE (10, 775) (HEAD (28, J) ,J=1 ,80) 

DO 130 1=34,36 

130 WRITE(I0,775)(HEAD(I,J),J=1,80) 

WRITE (10, 754)DTAUI,NISP,DTAUM,NMSP,DTAUF 
IF(ISH.EQ.l) GO TO 135 
WRITE (10, 756) 

IF(IRF.EQ.l) GO TO 135 
WRITE ( I 0 , 758 ) NFSP , ICL , NHSP 
135 WRITE (10,760) IFREQ , NTSP 

WRITE (10, 775) (HEAD(1,J), J=l,80) 

WRITE (10, 762) 

WRITE(I0,775)(HEAD(18,J),J=1,80) 

WRITE (10, 775 ) (HEAD(1 ,J) , J=l,80) 

INITIAL CONDITIONS 


MJ=M 

IRP=0 

INCL=1 

NIP=0 

JSH=1 

G2TIME=TG2/TREF 

G1TIME=TG1/TREF 

TMP=TMP/TREF 

X2TIME=TX2/TREF 

X1TIME=TX1/TREF 

TLEN1=TLEN1/12. 

TLEN2=TLEN2/12. 

TLEN3=TLEN3/12. 

DO 140 J=1,L 
140 EL(J)=EL(J)/12. 
TIN=TIN/TREF 
• T0N=T0N/3600. 
T0FF=T0FF/3600. 
DTAUI=DTAUI / 3600 . 
DTAUM=DTAUM/3600 . 
DTAUF=DTAUF/3600 , 
TAUI=TAUI / 3600 . 
TAUF=TAUF / 3 600 . 
TAFI=TAFI/3600. 
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DO 150 1=1, N 
• DO 150 J=1,M 
FDIST(I , J)=TIN 
IF(TIN.LT.TMP) GO TO 145 

H(I ,J)=CPL*TREF*FDIST(I , J)+DEL*(CPS*TMP*TREF/DES+HLAM) -TMP*TREF* 
1C PL 

145 H(I,J)=CPS*FDIST(I,J)*TREF 
150 CONTINUE 

INTERFACE NODE NUMBERING 

MM(1)=N0(1) 

DO 160 1=2, L 
MM(I)=MM(I-1)+N0(I)-1 
160 CONTINUE 

AKS=AK(L) 

INC=MM(L-1) 

DX=(TLEN2+TLEN3)/ (N*TLEN1) 

SET UP GAP NODE 

NODEH=IFIX (TLEN3/ (DX*TLEN1))+1 
N0DEG=N0DEH+1 
IF(TLEN3.EQ.0.) GO TO 165 
WRITE(I0,764)N0DEG 
GO TO 170 

165 WR1TE(I0, 764JN0DEH 

PRINT THE INITIAL TEMPERATURE PROFILE AND ENTHALPYS 

170 DO 175 1=1, N 
DO 175 J=1,M 
175 T(I,J)=FDIST(I ,J) 

DTAU=0. 

TAU=0 . 

WRITE (IO, 766) TAU 
DO 180 1=1, N 
DO 180 J=1,M 
180 TEO(I,J)=T(I,J)*TREF 

WRITE ( IO , 768 ) ((TEO(I,J) ,J=1,M) ,I=1,N) 

HSMP=CPS*TMP*TREF 

HLMP=DEL* (CPS'fTMP*TREF/DES+HLAM) 

HAFP= ( HSMP+HLMP ) / 2 . 

CALCULATION OF CONSTANTS 

LA=0 

DTAU=DTAUI 
200 IC0UNT=1 
MOW=MJ- 1 
NIP=NIP+1 


DETERMINE THE TIME STEP 
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IF(NIP .EQ .NISP+1 ) DTAU=DTAUM 
IF(JSH.EQ.3) GO TO 210 
IF(NIP.EQ.NMSP+1) DTAU=DTAUF 
210 IF(NIP , EQ . NHSP+1 ) DTAU=DTAUI 
LA=LA+1 
A1=1./DTAU 

IF(IRF.EQ.l) GO TO 215 
IF(NIP.EQ.NFSP) GO TO 687 

CONVERGENCE CRITERIA 

215 EPSI=0 .0001 

CALCULATION OF CONSTANTS 

NOW=N-l 
DO 220 K=1,L 

DY(K)=EL(K)/(TLEN1*(N0(K)-1)) 

C (K)=ALP (K ) / (2 . *TLEN 1*TLEN1*DX*DX ) 

D(K)=ALP(K)/(2.*TLEN1*TLEN1’ V DY(K)*DYCK)) 

A(K)=1./(A1+2.*(C(K)+D(K))) 

B(K)=A1-2.*(C(K)+D(K)) 

CPD(K)=AK(K)/ALP(K) 

220 CONTINUE 
NEW=L- 1 

IF(L.EQ.l) GO TO 230 
DO 225 K=1 ,NEW 

E(K)=1./(1./(D(K)*A(K))+(AK(K+1)*DY(K))/(AK(K)*DY(K+1)*A(K+1)* 

1D(K+1))) 

225 CONTINUE 
230 DO 240 1=1 ,N 
DO 240 J=1,MJ 
TO(I, J)=T(I , J) 

240 CONTINUE 

DO 250 1=1, N 
DO 250 J=INC,MJ 
HO(I,J)=H(I,J) 

250 CONTINUE 

TI.=TLEN1*TLEN1 

X1=1./(DX*DX) 

YL1=DY(L)*DY(L) 

YL2-DY f L" 1 )*DY (L- 1 ) 

E 1=2 . *TLEN 1*TLEN 1*D Y ( 1 ) *DY ( 1 ) / ( A ( 1 ) *ALP ( 1 ) ) 

Fl=4 . *H1*TLEN 1*DY ( 1 ) /AK ( 1 ) 

Gl=l./(El+Fl/2.) 

P 1=4 . *H2*TLEN1*DY (L) / AK (L) 

P2=l . /(I ./ (D(L)*A(L) )+Pl/2 . ) 

VAS=DTAU/(2.*TL*CPS) 

VAL=DTAU/ ( 2 . *TL*CPL) 

VA1=DTAU*AK (L- 1 ) / (2 . *CPS*TL*DY(L)*DY(L- 1 ) ) 

VA2=DTAU*AK ( L- 1 ) / ( 2 . *CPL*TL ,V DY ( L )*DY ( L - 1 ) ) 

VA3=DTAU*TREF/ (2 . *TL) 

VA4=AK(L-1)*DY(L-1)*X1/DY(L) 
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VA5=AK(L-1)/(DY(L-1)*DY(L)) 

VA9=2 . *H2*TLEN 1 /DY ( L) 

HEAT TERM INCLUSION 

IF(IHQ.EQ.1)Q1=Q 

IF( IHQ . EQ . 2 )Q1=QHEAT2 (TAU.DTAU ,TON ,TOFF ,QV) 

IF( IHQ . EQ . 3 )Q1=QHEAT3 (TAU , DTAU, TON, TOFF, QV) 

IF (IHQ . EQ . 4 )Q1=QHEAT4 (TAU.DTAU .TON ,TOFF ,QV) 

Q2=2. *3. 4121*144. *Q1 

CALCULATION OF CONSTANTS OF FINITE THICKNESS HEATER 

HE AT=Q2*DY ( I J ) *TLEN 1 *DY ( I J ) *TLEN 1 / ( TREF*AK (IJ)*EL(IJ)) 
FACT1=AK(IJ)*DY(IJ-1)/(AK(IJ-1)*DY(IJ)) 

FACT2=AK (I J+l )*DY( I J) / (AK( I J)*DY ( I J+l ) ) 

CALCULATION OF NEW TIME 

TAU=TAU+DTAU 

IF(L.NE.l) GO TO 295 
260 DO 270 1=1, N 
DO 270 J=2,MOW 
THOLD(I,J)=T(I,J) 

270 CONTINUE 

CONSTANT AMBIENT TEMPERATURE BOUNDARY CONDITIONS FOR 
ONE LAYER CASE ONLY 

DO 275 1=1 ,N 
T(I,1)=X1TIME 
T(I ,M)-X2TIME 
275 CONTINUE 

CALCULATION OF TEMPERATURES IN THE INTERIOR OF ONE LAYER 
CASE ONLY 

DO 285 J=2,MOW 

T(l,J)=A(l)*(C(l)*2.*(T(2,J)+TO(2,J))+D(l)*(T(l,J+l)+T(l,J-l)+ 

1T0(1,J+1)+T0(1,J-1))+B(1)*T0(1,J)) 

DO 280 1=2, NOW 

T(I,J)=A(1)*(C(1)*(T(I+1,J)+T(I-1,J)+T0(I+1,J)+T0(I“1,J))+ 

1D(1)*(T(I,J+1)+T(I,J-1)+T0(I,J+1)+T0(I,J-1))+B(1)*T0(I,J)) 

280 CONTINUE 

T(N,J)=A(l)*(C(l)*2.*(T(N-l,J)+TO(N-l,J))+D(l)*(T(N,J+l)+ 

1T(N,J-1)+T0(N,J+1)+T0(N,J-1))+B(1)*T0(N,J)) 

285 CONTINUE 


DUMAX=0 . 

CHECK CONVERGENCE FOR CONSTANT TEMPERATURE BOUNDARY CONDITIONS 
OF ONE LAYER CASE ONLY 
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DO 290 I»1 ,N 
DO 290 J=2,M0W 

DIF(I,J)=ABS(T(I,J)-THOLD(I,J)) 

290 DUMAX=AMAX1(DIF(I,J),DUMAX) 

GO TO 680 

295 IF((IBC1.EQ. 1).AND. (IBC2.EQ. 1 ) )G0 TO 315 

GAUSS -SEIDEL REITERATION 

300 DO 305 1*1, N 
DO 305 J=1,MJ 
TOLD (I ,J)=T(I ,J) 

305 CONTINUE 

CALCULATION OF TEMPERATURES AT THE INNER-AMBIENT INTERFACE 

T(l, 1)=G1*(DY(1)*DY(1)*2 .*(T(2, l)+TO(2i l))*Xl+2 .*(T(1 ,2) *■ 
1T0(1 > 2))+(B(1)/D(1)-F1/2.)*T0(1,1)+F1*G1TIME) 

DO 310 1=2, NOW 

T(I , 1)=G1*(DY(1)*DY(!)*(T(I+1 ,1)+T(I-1 ,1)+T0(I+1 ,1)+T0(I-1,1)) 
1*X1+2.*(T(I ,2)+T0(I ,2) )+(B(l)/D(l)-Fl/2..)*T0(I , 1)+F1*G1TIME) 
310 CONTINUE 

T(N,1)=G1*(DY(1)*DY(1)*2>(T(N-1,1)4T0(N-1,1))*X1+2.*(T(N,2)+ 

1T0(N,2))+(B(1)/D(1)-F1/2.)*T0(N,1)+F1*G1TIME) 

GO TO 330 

315 DO 320 1=1, N 
DO 320 J=2 , MOW 
PHOLD (I,J)=T(I,J) 

320 CONTINUE 

CONVECTION BOUNDARY CONDITIONS FOR ONE LAYER CASE ONLY 

DO 325 1=1 ,N 
T(I,1)=X1TIME 
T(I,M)=X2TIME 
325 CONTINUE 

330 LIST=2 


CALCULATION OF TEMPERATURES IN THE INTERIOR OF THE LAYER 
AND AT THE INTERFACE BETWEEN LAYERS 


DO 600 K=1,L 
NN=MM(K) 

NNEW=NN-1 

INC1=INC+1 

IF(IH.EQ. 3) GO TO 335 
GO TO 340 

335 IFCNN.EQ.N02) GO TO 445 
340 IF(LIST.GT.NNEW) GO TO 360 
DO 355 J=LIST,NNEW 
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IFUG.EQ.l) GO TO 345 
IF(J.EQ.INCl) GO TO 530 

345 T(1 ,J)»A(K)*(C(K)*2 .*(T(2 ,J)+T0(2 ,J) )+D(K)*(T(l ,J+1)+ 
1T(1 I J-1)+TO(1,J+1)+TO(1,J-1))+B(K)*TO(1 > J)) 

DO 350 1*2, NOW 

T(I,J)«A(K)*(C(K)*(T(I+1,J)+T(I-1, J)+TO(I+l,J)+TQ(I-l,J))+ 
1D(K)*(T(I,J+1)+T(I,J-1)+TO(I,J+1)+TO(I,J-1))+B(K)*TO(I,J)) 

350 CONTINUE 

T(N,J)=A(K)*(C(K)*2.*(T(N-l,J)+TOCN-l,J))+D(K)*(T(N,J+l)+ 

1T(N,J-1)+TO(N,J+1)+TO(N,J-1))+B(K)*TO(N,J)) 

355 CONTINUE 
360 LIST=NN+1 

IF(IG.EQ.l) GO TO 375 

IFCISH.EQ.1) GO TO 370 

IF(IRF.EQ.l) GO TO 370 

IF(IRP.EQ.l) GO TO 365 

IF(NIP.GE.NFSP) GO TO 370 

IF( (JSH.EQ.2) .AND. (NN.EQ . INC)) GO TO 610 

GO TO 370 

365 IF((JSH.EQ.3) .AND. (NN.EQ. INC)) GO TO 610 
370 IF(NN.EQ.INC) GO TO 475 
375 IF(NN.EQ.MJ) GO TO 605 
IF(IH.EQ.3) GO TO 380 
IF(IH,EQ.2) GO TO 385 
IF(IH.EQ.l) GO TO 395 
380 IF(NN.EQ.NOl) GO TO 430 
GO TO 395 

385 IF (NN.EQ. NODE) GO TO 390 
GO TO 395 

390 IF(TLEN3.EQ.O.) GO TO 415 

395 T(1 ,NN)=E(K)*( (DY(K)*DY(K)/ (DX*DX)+DY(K+1)*DY(K)*AK(K+1)/ 

1(DX*DX*AK(K)))*2.*(T(2,NN)+T0(2,NN))+2.*(T(1,NN“1)+T0(1,NN-1))+ 

22.*AK(K+l)*DY(K)*(T(l,NN+l)+TO(l,NN+l))/(AK(K)*DY(K+l))+(B(K)/ 

3D(K)+AK(K+l)*DY(K)' fr B(K+l)/(AK(K)*D(K+l)*DY(K+l)))*T0(l,NN)) 

DO 410 1=2, NOW 
IFCIH.EQ.2) GO TO 400 
IF(IH.EQ.l) GO TO 405 
IF(IH.EQ.3) GO TO 405 

400 IF((NN. EQ. NODE). AND. (I.EQ.NODEG)) GO TO 420 

405 T(I,NN)=E(K)*((DY(K)*DY(K)/(DX*DX)+DY(K+1)*DY(K)*AK(K+1)/ 

l(DX*DX*AK(K)))*(T(I+l,NN)+T(I-l,NN)+T0(I+l,NN)+T0(I-l,NN))+2.* 

2(T(I,NN*l)+TO(I,NN-l))+2.*AK(K+l)*DY(K)*V(T(I,NN+l)+TO(I,NN+l))/ 

3(AK(K)*DY(K+1))+(B(K)/D(K)+AK(K+1)*DY(K)*B(K+1)/(AK(K)* 

4D (K+l )*DY (K+l ) ) )*T0 ( I , NN) ) 

410 CONTINUE 

T(N,NN)=E(K)*((DY(K)*DY(K)/(DX*DX)+DY(K+1)*DY(K)*AK(K+1)/ 
1(DX*DX*AK(K) ) )*2 .*(T(N-1 ,NN)+T0(N-1 ,NN) )+2 .*(T(N,NN-1)+ 
2T0(N,NN-1))+2.*AK(K+1)*DY(K)*(T(N,NN+1)+T0(N,NN+1))/(AK(K)* 
3DY(K+1) )+(B(K)/D(K)+AK(K+l)*DY(K)*B(K+l)/(AK(K)*D(K+l)* 

4DY (K+l ) ) )*T0 (N , NN) ) 

GO TO 600 

POINT HEATER 
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415 T(1,N0DE)»E(L1)*((DY(L1)*DY(L1)/(DX*DX)+DY(L2)*DY(L1)*AK(L2)/ 
1(DX*DX*AK(L1)))*2.*(T(2,N0DE)+T0(2,N0DE))+2.*(T(1,N0D£-1)+ 

2T0 ( i , NODE - 1 ) ) +2 . *AK ( L2 ) *DY ( L 1 ) * (T ( 1 , N0DE+ 1 ) +TO ( 1 , N0DE+ 1 ) ) / 
3(AK(L1)*DY(L2)H(B(L1)/D(L1)+AK(L2)*DY(L1)*B(L2)/(AK(L1)* 
4D(L2)*DY(L2)))nO(l,NODE)+4.nLENl*DY(Ll)*Ql*3. 4121*144,/ 
5(TREF*AK(L1))) 

420 DO 425 I=NODEG,NOW 

T(I,N0DE)*E(L1)*((DY(L1)*DY(L1)/(DX*DX)+DY(L2)*DY(L1)*AK(L2)/ 

1 (DX*DX*AK(L1 ) ) )* (T ( 1+1 , NODE ) +T ( I - 1 , NODE ) +TO (1+1 .NODE )+ 
2T0(I-1,N0DE))+2.*(T(I,N0DE-1)+T0(I,N0DE-1))+2.*AK(L2)*DY(L1)* 
3(T(I,N0DE+1)+T0(I,N0DE+1))/(AK(L1)*DY(L2))+(B(L1)/D(L1)+ 
4AK(L2)*DY(Ll)*B(L2)/(AK(Ll)*D(L2)*DY(L2)))*T0(I,N0DE)+4.* 
STLEN1*DY (LI ) *Q1*3 . 4 12 1*144 . / (TREF*AK ( LI ) ) ) 

425 CONTINUE 

T(N,N0DE)=>E(L1)*((DY(L1)*DY(L1)/(DX*DX)+DY(L2)*DY(L1)*AK(L2)/ 

1 (DX*DX*AK (LI ) ) )*2 . * (T(N- 1 , NODE )+TO (N - 1 , NODE ) ) +2 . * (T ( N , NODE- 1 ) + 
2TO(N, NODE-1) )+2.*AK(L2)*DY(Ll)*(T(N,NODE+l)+TO(N,NODE+l))/ 
3(AK(L1)*DY(L2) )+(B(L,l)/D(Ll)+AK(L2)*DY(Ll)*B(L2)/(AK(Ll)* 
4D(L2)*DY(L2)))*T0(N,N0DE)4*.*TLEN1*DY(L1)*Q1*3. 4121*144./ 
5(TREF*AK(L1)) ) 

GO TO 600 

FINITE THICKNESS HEATER 

430 IF(TLEN3.NE.O. ) Q7=0. 

IFCTLEN3.EQ.O.) Q7=HEAT 

T(l,N01)oE(IJ-l)*((DY(IJ-i)*DY(IJ-l)/(DX*DX)+DY(IJ)*DY(IJ)* 

1FACT1/(DX*DX))*2.*(T(2,N01)+T0(2,N01))+2.*(T(1.N01-1)+ 

2T0(1,N01-1))+2.*FACT1*(T(1,N01+1)+T0(1,N01+1))+(B(IJ-1)/D(IJ-1)+ 

3B(IJ)*FACT1/D(IJ))*T0(1,N01)+Q7*FACT1) 

DO 435 I«2,N0W 
IF(I.LT.NODEG) Q8=0. 

IF(I.GE.NODEG) Q8=HEAT 

T(I,N01)=E(IJ-1)*((DY(IJ-1)*DY(IJ-1)/(DX*DX)+DY(IJ)*DY(IJ)* 
1FACT1/(DX*DX))*(T(I+1,N01)+T(l-1,N01)+T0(l+1,N01)+T0(l-1,N01))+ 
22.*(T(I ,NQ1"1)+T0(I .NOl-l) )+2 .*FACT1*(T(I ,N01+1)+T0(I ,N01+1))+ 
3(B(IJ-1)/D(IJ-1)+B(IJ)*FACT1/D(IJ))*T0(I,N01)+Q8*FACT1) 

435 CONTINUE 

T(N,N01)=E(IJ-1)*((DY(IJ-1)*DY(IJ-1)/(DX*DX)+DY(IJ)*DY(IJ)* 

1FACT1/(DX*DX))*2.*(T(N-1,N01)+T0(N-1,N01))+2.*(T(N,N01-1)+ 

2T0(N,N01-l))+2.*FACTl*(T(N,NOl+l)+T0(N,N01+l))+(B(IJ-l)/D(IJ-l)+ 

3B(IJ)*FACT1/D(IJ))*T0(N,N01)+HEAT*FACT1) 

GO TO 600 

445 IF ( (N02 -NO 1 ) . EQ . 1 ) GO TO 465 

MORE THAN TWO NODES IN THE HEATER 

450 NFW=N01+1 
NNFW=N02-1 
DO 460 IL=?^y.NNFW 
IF(TLEN3.£0, Q3=Q2 
IF(TLEN3.NE.O. ) Q3=0. 

T(1,IL)=A(IJ)*(C(IJ)*2.*(T(2,IL)+T0(2,IL))+D(IJ)*(T(1,IL+1)+ 
1T(1»IL-1)+T0(1 , IL+1)+T0(1 , IL-1) )+B(IJ)*TO(l , IL)+Q3/ (EL(IJ)*TREF* 
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22.*CPD(IJ))) 

DO 455 I«2,N0W 
IF(I.LT.NODEG) Q4*0. 

IF(I.GE.NODEG) Q4«Q2 

T(I ,IL)"A(IJ)*(C(IJ)*(T(I+1 ,IL)+T(I-1,IL)+T0(I+X, IL)+T0(I*1 ,IL))+ 
lD(IJ)*(Ta,IL+l)+Ta,IL-l)+T0(I,IL4-l)+T0(I,IL-l))+B(IJ)*T0(I,IL)+ 
2Q4/ (EL(I J)*TREF*2 .*CPD(I J) ) ) 

455 CONTINUE 

T(N, IL)-A( I J)*(C(IJ)*2 >(T(N-1 , IL)+TO(N-l , IL) )+D(IJ)*(T(N, IL+1J+ 
1T(N , IL- 1 )+TO (N , I L+l )+TO(N , I L- 1 ) )+B C I J)*TO(N , IL)+Q2/ (TREF*EL( I J )* 
22.*CPD(IJ))) 

460 CONTINUE 

465 IF(TLEN3.EQ.O.) Q5-HEAT 
IF(TLEN3.NE .0. ) Q5»0. 

T(1 I N02)«E(IJ)*((DY(IJ)*DY(IJ)/(DX*DX)+DY(IJ+1)*DY<IJ+1)*FACT2/ 
l(DX*DX))*2.*(T(2,N02)+T0(2,NO2))+2.*(T(l,NO2-l)+T0(l,N02-l))+2.* 
2FACT2*(T(l,NO2+l)+TO(l,N02+l))+(B(IJ)/D(IJ)+B(IJ+l)*FACT2/ 
3D(IJ+1))*T0(1,N02)+Q5) 

DO 470 I®2,N0W 
IF (I .LT.N0DE6) Q6*0. 

IF ( I .GE.NODEG) Q6»HEAT 

T(I,N02)=E(IJ)*((DY(IJ)*DY(IJ)/(DX*DX)+DY(IJ+1)*DY(IJ+1)*FACT2/ 
l(DX*DX))*(T(I+l,NO2)+T(I-l,NO2)+T0(I+l,NO2)+TO(I-l,NO2))+2.* 
2(T(I,N02«l)+TO(I,N02-l))+2>FACT2*(T(I,N02+l)+T0(I,N02+l))+(BCIJ) 
3/D(IJ)+B(I J+1)*FACT2/D( I J+l) )*TO( I f N02)+Q6) 

470 CONTINUE 

T (N , N02 )=E ( I J.)* ( (DY ( I J)*DY ( I J) / (DX*DX)+DY ( I J+l )*DY ( I J+ 1 )*FACT2/ 
l(DX*DX))*2.*(T(N-l,N02)+T0(N-l,NO2))+2.*(T(N,N02-l)+TO(N,N02-l))+ 
22.*FACT2*(T(N,N02+1)+T0(N,N02+1) )+(B(IJ)/D(IJ)+B(IJ+l)*FACT2/ 

3D ( I J+l ) ) *T0 (N , N02 )+HEAT) 

LIST=NN+1 
GO TO 600 

CALCULATION OF ENTHALPYS AND TEMPERATURES AT THE SHIELD- 
WATER INTERFACE 

475 IF(ICOUNT.NE.l) GO TO 480 

CALL CHANGE (HO (2, INC) t HOC 1, INC) ,HO(2, INC) ,HO( 1 , INC1) ,HO(l , INC) , 
1IP0 ,AK01 ,AK02 ,AK03 ,AK04) 

CONST ( 1 , INC)=AK02*X1*2 . * (TO (2 , INC ) -TO(l , INC ) )+2 . *AK03* (TO ( 1 , INC+1 ) 
1-T0(1,INC))/YL1 

480 CALL CHANGE (H( 2, INC) ,H(1, INC) ,H(2, INC) , H ( 1 , INC 1 ) ,H(1 ,INC) , IP, 
1AK1,AK2,AK3,AK4) 

VA6=(TMP*TREF- (DEL* (CPS*TMP*TREF/DES+HLAM) /CPL) ) /TREF 
VK1=(AK1+AK2)*X1 
VK2= (AK3+AK4 ) / YL1 
VK3=2 . *AK3/YL1 
VK4=2 .*AK4/YL1 
GO TO(483,485 ,487) , IP 
C ICE 

^ 483 H(1 J INC)=1./(VA1/(D(L-1)*A(L-1))+VAS*(VK1+VK3)+1.)*(H0(1,INC)+ 
1VA3* ( (T(2 , INC)+TO (2 , INC ) )*2 . *VA4+X1*2 . *AK2*T(2, INC )+2 . *VA5* ( 

2T(1 , INC-1 )+TO( 1 , INC-1) )+VK3*T( 1 , INC+l)+CONST( 1 , INC)+B (L-1)*VA5* 
3TO(l,INC)/D(L-l))) 
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T( 1 , INC)«K( l , INC)/ (CPS*TREF) 

IF(H(l f INC) . LE.HSMP) GO TO 490 
IP*1 

C MELTING POINT 

485 H( 1 , INC)*HO( 1 , INC)+VA3*(- (VA5/ (D(L- 1)*A(L-1) )+VKl+VK3)*TMP+( 

1T(2 ,INC)+TO(2 , INC) )*2 .*VA4+X1*AK2*2 .*T(2 ,INC)+2 .*VAS*(T(1 ,INC-1)+ 
2TO( 1 , INC- 1 ) )+VK3*T( 1 , INC+1 )+CONST( 1 , INC )+B(L- 1 )*VA5*TO( 1 , INC)/ 
3D(L-1)) 

T(1,INC)«TMP 

IF(H(1, INC). GT.HSMP.AND.HU, INC). LT.HLMP) GO TO 490 
IF(IP.EQ. 1.0R.IP.EQ.3) GO TO 490 
C LIQUID WATER 

487 H(1 , INC)® I. / (VA2/ (DCL-1 )*A(L-1) )+VAL*(VKl+VK3)+l . )*(HO( 1 , INC)+ 
1VA3*< (T(2 , INC)+TO(2 , INC) )*2 .*VA4+X1*2 .*AK2*T(2 , INC)+2 .*VA5*( 

2T( 1 , INC- 1 )+TO (1 , INC- 1 ) )+VK3*T( 1 , INC+1 )+C0NST( l , INC) - (VA5/ (D( L- 1)* 
3A(L- 1 ) )+VK 1+VK3 ) *VA6+B ( L* 1 )*VA5*TO (1 . INC ) /D(L- 1 ) ) ) 

T(1 , INC)«H(1 , INC)/(TREF*CPL) -DEL*(CPS*TMP*TREF/DES+HLAM)/ (CPL*TREF 
1)+TMP 

IF(H(1,INC).GE.HLMP) GO TO 490 
IP® 3 

GO TO 483 

490 DO 510 I m 2,N0W 

IF(ICOUNT.NE.l) GO TO 495 

CALL CHANGE (HO ( I - 1 , ItyG ) ,HO(I .INC) ,H0(I+1 ,INC) ,HO(I ,INC1) , 

1H0( I , INC ) , IPO , AKO 1 , AK02 ,AKC3, AK04 ) 

C0NST( I , INC)=AK01*Xi* (T0( I - 1 , INC) -T0( I , INC) ) -AK02*X1* (TO(I , INC ) • 
1T0 ( 1+ 1 , INC ) ) +2 . *AK03* (TO ( I , I NC+ 1 ) -TO ( I , INC ) ) / YL1 
495 CALL CHANGE (H( I -1, INC) ,H(I , INC) ,H(I+1,INC) ,H(I ,INC1) ,H(I ,INC) , IF' 
1AK 1 , AK2 , AK3 , AK4 ) 

VK1=(AK1+AK2)*X1 
VK2 E (AK3+AK4)/YL1 
VK3=2 . *AK3/YL1 
VK4=2.*AK4/YL1 
GO TO(503 ,505 ,507 ) , IP 
C ICE 

503 H(I,INC)=1./(VA1/(D(L-1)*A(L-1))+VAS*(VK1+VK3)+1.)*(H0(I,INC)+ 
1VA3*( (T( 1+1 , INC )+T( I - 1 , INC)+TO( 1+1 , INC)+TO( I - 1 , INC) )*VA4+X1* (AK2* 
2T(I+1,INC)+AK1»W-1, INC) )+2.*VA5*(T(I,INC-l)+TO(I, INC-1) )+VK3* 

3T ( I , INC+ 1 ) +CONST (I , INC )+B ( L- 1 )*VA5*T0 (I,INC)/D(L-1))) 

T(I, INC)=H(I , INC)/ (CPS*TREF) 

IF(H(I,INC) .LE.HSMP) GO TO 510 
IP=1 

C MELTING POINT 

505 H(I,INC)«=H0(I,INC)+VA3*(-(VA5/(D(L-l)*A(L-l))+VKl+VK3)nMP+( 

1T( 1+1 , INC)+T (I - 1 , INC )+T0 ( 1+1 , INC )+T0 ( I - 1 , INC ) )*VA4+X1*(AK2* 

2T( 1+1 , INC)+AK1*T(I - 1 , INC) )+2 . *VA5* (T( I , INC- 1 )+T0 ( I , INC- 1) )+VK3* 
3T(I,INC+1)+C0NST(I,INC)+B(L-1)*VA5*T0(I,INC)/D(L-1)) 

T(I,ING)*TMP 

IF(H(I,II-: GT.HSMP.AND.HU, INC). LT.HLMP) GO TO 510 

IFCIP.EQ- ’ .0R.IP.EQ.3) GO TO 510 
C LIQUID VAKR 

507 H(I,INC)=1./(VA2/(D(L-1)*A(L-1))+VAL*(VK1+VK3)+1.)*(H0(I,INC)+VA3* 
1 ( (T ( 1+1 , INC ) +T ( I - 1 , INC )+T0 ( 1+ 1 , INC ) +T0 ( I - 1 , INC ) ) *VA4+X1* ( AK2* 
2T(I+1,INC)+AK1*T<I-1,INC))+2.*VA5*(T(I,INC-1)+T0(I,INC-1))+VK3* 
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3T(I,INC+1)+C0NST(I,INC)-(VA5/(D(L-1)*A(L-1))+VK1+VK3)*VA6+B(L-1)* 

4VA5*TO(I,INC)/D(L-l))) 

T(I,INC)-H(I,INC)/<TREF*CPL)-DEL*(CPS*TMP*TREF/DES+HLAM)/(CPL*TREF 

1)+TMP 

IF (H( I , INC) .GE.HLMP) GO TO 510 
IP**3 

GO TO 503 
510 CONTINUE 

IFCICOUNT.NE.I) GO TO 520 

CALL CHA.NGE(HQ(N-1 ,INC) ,KO(N,INC) ,H0(N-1,INC) ,H0(N,INC1) ,HO(N,INC) 
1 , IPO , AK01 , AK02 i AK03 , AK04 ) 

CONST(N . INC ) «AK0 J *X 1*2 * (TO (N- 1 , INC ) -TO (N , INC ) )+2 . *AK03* ( 

1T0 (N , INC+ 1 ) -TO(N , INC ) ) / YL1 

520 CALL CHANGE (H(N-1, INC) ,H(N, INC) ,H(N- 1 , INC) ,H(N , INC1) ,H(N, INC) , IP, 
2AK 1 , AK2 , AK3 , AK4 ) 

VK1=(AK1+AK2)*X1 
VK2 = (AK3+AK4)/YL1 
VK3“2.*AK3/YL1 
VK4=2.*AK4/YL1 
GO TO(523 ,525 ,527) ,IP 
C ICE 

523 H(N, INC)*1 . / (VAX/ (D(L-1 )*A(L- 1 ) )+VAS*(VKl+VK3)+l . )*(HO{N, INC)+ 
).VA3*( (T(N-1 , INC)+T0(N-1 ,INC) )*2 .*VA4+X1*2 ,*AK1*T(N-1 ,INC)+2 .*VA5* 
2(T(N,INC-1)+T0(I,INC-1))+VK3*T(N,INC+1)+C0NST(N,INC)+B(L-1)*VAS* 
4T0(N,INC)/D(L-1))) 

T(N,INC)«H(N,INC)/(TREF*CPS) 

IF(H(N,INC) .LE.HSMP) GO TO 600 
IP-=1 

C MELTING POINT 

525 H(N,INC)=H0(N,INC)+VA3*(-(VA5/(D(L-l)*A(L-l))+VKl+VK3)*TMP+< 

1T(N-1 ,INC)+T0(N* 1 , INC) )*2 .*VA4+X1*AK1*2 .*T(N-1 , INC)+2 .*VA5*( 
2T(N,INC-l)+T0(N,INC-l))+VK3*T(N,INC-fl)+C0NST(N,INC)+B(L-l)*VA5* 
3T0(N,INC)/D(L-1) ) 

T(N,INC)=*TMP 

IF(H(N, INC) .GT>HSMP. AND.H(N, INC) .LT.HLMP) GO TO 600 
IFCIP.EQ . 1 .0R.IP.EQ.3) GO TO 600 
C LIQUID WATER 

527 H(N,INC)=1./(VA2/(D(L-1)*A(L-1))+VAL*(VK1+VK3)+1.)*(H0(N,INC)+VA3* 
1 ( (T(N-1 , INC)+TO(N- 1 , INC) )*2 . *VA4+X1*2 .*AK1*T(N-1 , INC)+2 .*VA5*( 

2T(N , INC- 1 )+T0 (N , INC- 1 ) )+VK3*T(N , INC+1 )+CONST(N , INC) - (VA5/ (D(L- 1 )* 
3A(L- 1 ) )+VKl+VK3 )*VA6+B (L- 1 )*VA5*T0 (N , INC ) /D (L- 1 ) ) ) 

T (N , INC )=H (N , INC ) / (TREF*CPL) -DEL* (CPS*TMP*TR£F/DES+HLAM) / (CPL*TREF 
D+TMP 

IF(H(N, INC). GE.HLMP) GO TO 600 
IP=3 

GO TO 523 

CALCULATION OF ENTHALPYS AND TEMPERATURES IN THE INTERIOR OF 
THE ICE-WATER LAYER 

530 DO 590 J=INC1,M0W 

IF(ICOUNT.NE.l) GO TO 535 

CALL CHANGE(H0(2, J) ,H0(1, J) ,H0(2,J) ,H0(1,J+1) ,H0(1,J-1) ,IPO, 

1 AKO 1 , AK02 , AKO 3 , AK04 ) 
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CONST(l,J)=AKO2<' r Xl*2.*(TO(2,J[)-TO(l,J))+(AK03*(T0(l l J+l)-TO(l,J) 1 
1-AK04*(TO(1,J)-TO(1,J-1)))/YL1 

535 CALL CHANGE (H( 2, J) ,H(1 , J) ,H(2, J) ,H(1 , J+l) ,H(1, J-l) ,IP,AK1 ,AK2, 
1AK3.AK4) 

VK 1= ( AK 1+AK2 ) *X 1 

VK2= ( AK3+AK4 ) / YL1 

VK3=2.*AK3/YL1 

VK4=2.*AK4/YL1 

GO TO (543, 545, 54 7), IP 

ICE 

543 H(1 , J) e l. /(I .+VAS*(VKl+VK2))*(HO(l, J)+VA3*(Xl rt 2 .^AK2*T(2, J)+( 
1AK3*T( 1 , J+l )+AK4*T( 1 , J- 1 ) )/YLl+CONST( 1 , J) ) ) 

T(1,J)=H(1,J)/(CPS*TREF) 

IF(H(1,J) .LE.HSMP) GO TO 550 
IP=1 

MELTING POINT 

545 H ( 1 , J ) =HQ ( 1 , J ) +VA3* ( - ( VK 1 +VK2 )*TMP+X 1*2 . *AK2*T ( 2 , J )+ ( AK3*T ( 1 , J+ 1 ) + 
1 AK4*T (1,J-1))/ YL 1+CONST ( 1 , J ) ) 

T(1,J)=TMP 

IF(H(1,J) ,GT.HSMP.AND.H(1,J) .LT.HLMP) GO TO 550 
IF(IP.EQ. l.OR. IP.EQ.3) GO TO 550 
LIQUID WATER 

547 H(l,J)=l./(l. +VAL*(VK1+VK2 ) )* (HO ( 1 , J )+VA3* ( - ( VK1+VK2)*VA6+X1* 
12.*AK2*T(2,J)+(AK3*T(1,J+1)+AK4*T(1,J-1))/YL1+C0NST(1,J))) 

T( 1 , J)=H( 1 , J) / (CPL*TREF) -DEL*(CPS*TMP*TREF/DES+HLAM) / (CPL*TREF)+ 
1TMP 

IF(H(1,J) .GE.HLMP) GO TO 550 
IP=3 

GO TO 543 

550 DO 570 1=2, NOW 

IF(ICOUNT.NE.l) GO TO 555 

CALL CHANGE (H0( I -1,J) ,HO(I ,J) ,H0(I+1,J) ,HO(I , J+l) ,HO(I , J-l) , 

1 I PO , AKO 1 , AK02 , AK03 , AK04 ) 

CONST(I,J)=AK01*X1*(TO(I-1,J)-TO(I,J))-AK02*X1*(TO(I,J)-TO(I+1,J)) 
1+(AK03*(T0(I,J+1)-T0(I,J))-AK04*(T0(I,J)-T0(I,J-1)))/YL1 
555 CALL CHANGE (H(!-1,J) ,H(I , J) ,H(I+1 , J) ,H(I, J+l) ,H(I , J-l) ,IP , 

1AK1 , AK2 , AK3 , AK4 ) 

VK 1= ( AK 1+AK2 ) *X i 

VK2=(AK3+AK4)/YL1 

VK3=2.*AK3/YL1 

VK4=2.*AK4/YL1 

GO T0(563,565 ,567) ,IP 

ICE 

563 H(I,J)=l./(l.+VAS*(VHi+VK2))*(H0(l,J)+VA3*(Xl*(AK2*T(I+l,J)+AKl* 
1TT-1,J))+(AK3*T(I,J+1)+AK4*T(I,J-1))/YL1+C0NST(I,J))) 
T(I,J)=H(I,J)/(CPS*TREF) 

IF(H(I,J). LE.HSMP) GO TO 570 
IP=1 

MELTING POINT 

565 H(I,J)=HO(I, J)+VA3*(-(VK1+VK2)*TMP+X1*(AK2*T(I+1,J)+AK1*T(I-1,J))+ 
1 ( AK3*T ( I , J+l )+AK4*T ( I , J- 1 ) ) / YL1+C0NST ( I , J ) ) 

T(I , J)=TMP 

IF(H(I,J).GT.HSMP. AND. H(I,J). LT.HLMP) GO TO 570 
IF(IP.EQ. l.OR. IP.EQ.3) GO TO 570 
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C LIQUID WATER 

567 H(I , J)=l ./(l .+VAL*(VK1+VK2) )*(HO(I , J)+VA3*(-(VKl+m)*VA6+Xl*(AK2* 
1T(I+1, J)+AK1*T(I-1,J))+(AK3*T(I,J+1)+AK4*T(I , J-1))/YL1+C0NST(I , J)) 
2) 

TH , J)=H(I ,J)/(TREF*CPL)-DEL*(CPS*TMP*TREF/DES+HLAM)/(CPL*TREF)+ 
1TMP 

IF(H(I,J) . GE.HLMP) GO TO 570 
IP=3 

GO TO 563 
570 CONTINUE 

IF(ICQUNT.NE.l) GO TO 575 

CALL CHANGE (HO(N- 1 ,J) ,HQ(N, J) ,H0(N-1 , J) ,HO(N,J+l) ,H0(N, J-l) , 

1 IPO , AK01 , AK02 , AK03 , AK04 ) 

CONST(N,J)=AK01*X1*2.*(TO(N-1,J)-TO(N,J))+(AK03*(TO(N,J+1)-TO(N,J) 

1)-AK04*(T0(N,J)-T0(N,J-1)))/YL1 

575 CALL CHANGE(H(N-1 , J) ,H(N, J) ,H(N-1, J) ,H<N,J+1) ,H(N, J-l) ,IP,AK1, 
1AK2,AK3,AK4) 

VK1=(AK1+AK2)*X1 
VK2= ( AK3+AK4 ) / YL1 
VK3=2.*AK3/YL1 
VK4=2.*AK4/YL1 
GO TO (583, 585, 587) ,IP 
C ICE 

583 H(N,J)=1./(1.+VAS*(VK1+VK2))*(H0(N,J)+VA3*(X1*AK1*2.*T(N-1,J)+ 

1 (AK3*T (N, J+l )+AK4*T(N , J- 1 ) ) /YL1+CONST (N , J) ) ) 

T(N,J)=H(N, J)/(CPS*TREF) 

IF(H(N, J) .LE.HSMP) GO TO 590 
IP=1 

C MELTING POINT 

585 H(N,J)=H0(N,J)+VA3*(-(VK1+VK2)*TMP+X1*2.*AK1*T(N-1,J)+(AK3* 

1T(N , J+l )+AK4*T(N , J- 1 ) ) /YL1+C0NST(N, J) ) 

T(N,J)=TMP 

IF(H(N,J) .GT.HSMP.AND.H(N,J) .LT.HLMP) GO TO 590 
IF(IP.EQ. 1 .OR. IP.EQ.3) $0 TO 590 
C LIQUID WATER 

587 H(N, J)-l . /(l.+VAL*(VKl+VK2) )*(HO(N, J)+VA3*(-(VK1+VK2)*VA6+X1*2 .* 
1AK1*T(N- 1 , J)+(AK3*T(N , J+l )+AK4*T(N , J- 1 ) ) /YL1+C0NST(N , J) ) ) 

T(N, J)”+(N, J)/ (TREF*CPL) -DEL* (CPS*TMP*TREF/DES+HLAM) / (CPL*TREF)+ 
1TMP 

IF(H(N,J) .GE.HLMP) GO TO 590 
IP=3 

GO TO 583 
590 CONTINUE 
600 CONTINUE 

605 IF((IBC1.EQ. 1) .AND. (IBC2.EQ.1))G0 TO 672 
IF(IG.EQ. 2) GO TO 615 
C 

C CALCULATION OF TEMPERATURES AT THE OUTER-AMBIENT INTERFACE 
C WITHOUT PHASE CHANGE 
C 

610 T(1,MJ)=P2*(YL1*2.*(T(2,MJ)+T0(2,MJ))*X1+2.*(T(1,MJ-1)+ 
1T0(1,MJ-1))+(B(L)/D(L)-Pi/2.)*T0(1,MJ)+P1*G2TIME) 

DO 612 1=2, NOW 

T(I,MJ)=P2*(YL1*(T(I+1,MJ)+T(I-1,MJ)+T0(I+1,MJ)+T0(I-1,MJ))*X1 
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1+2 .*(T(I ,MJ-l)+TO(I ,MJ-l))+(B(L)/D(L)*Pl/2 . )*TO(I ,MJ)+P1*G2TIME) 
612 CONTINUE 

T(N,MJ)=P2*(YLl*2.*(T(N-l,MJ)+TO(N-l,MJ))*Xl+2. 

l*(T(N,MJ-l)+T0(N,MJ-l))+(Ba)/D(L)-Pl/2.)*T0(N,MJ)+Pl*G2TIME) 

GO TO 677 

CALCULATION OF TEMPERATURES AND ENTHALPYS AT THE OUTER- 
AMBIENT INTERFACE 

615 IF(ICOUNT.NE.I) GO TO 617 

CALL CHANGE (H0(2,M) ,HO(l ,M) ,HO(2,M) ,HO(l ,M) ,HO(l ,M-1) ,IPO, 

1 AKO 1 , AK02 , AK03 , AK04 ) 

CONST ( 1 , M)=-AK02*X1* V 2 . * (TO ( 2 , M) -TO ( 1 , M) )+2 . *AK04* (TO( 1 , M) - 
1T0(1 ,M-1))/YL1 

617 CALL CHANGE (H (2 ,M) ,H(1 ,M) ,H(2 ,M) ,H(1 ,M) ,H(1 ,M-1) , IP ,AK1 , 
1AK2.AK3.AK4) 

VK1=(AK1+AK2)*X1 

VK2=(AK3+AK4)/YL1 

VK3=2,*AK3/YL1 

VK4=2.*AK4/YL1 

GO TO(623, 625,627), IP 

ICE 

623 H(1,M)=1./(1.+VAS*(VK1+VK4+VA9))*(H0(1,M)-VA3^(-X1*2.*AK2*T(2,M)- 
1VK4*T(1 ,M-l)+CONST(l ,M) -VA9*(G?TIME*2 . -TOC1 ,M) ) ) ) 
T(1,M)=H(1,M)/(CPS*TREF) 

IF(H(1,M).LE.HSMP) GO TO 630 
IP=1 

MELTING POINT 

625 H(1 ,M)=HO(l ,M) -VA3*( (VK1+VK4+VA9)*TMP-X1*2 .’'AK2*T(2 ,M) -VK4* 

IT ( 1 , M - 1 ) +CONST ( 1 , M ) -VA9* (G2TIME*2 . -TO ( 1 , M) ) ) 

T(1 ,M)=TMP 

IF(H(1,M) .GT.HSMP.AND.H(1,M) .LT.HLMP) GO TO 630 
IF(IP.EQ. 1 .OR. IP.EQ. 3) GO TO 630 
LIQUID WATER 

627 H(1 ,M)=1. / (1 . +VAL* (VK1+VK4+VA9 ) ) * (HO ( 1 ,M) -VA3* ( (VK1+VK4+VA9 
lVA6-Xl ,lf AK2*2 . *T(2 ,M) -VK4*T( 1 ,M- 1 )+CONST( 1 , M) -VA9*(G2TIME*2 . - 
1T0(1,M)))) 

T(1 ,M)=H(1 ,M)/ (TREF*CPL) -DEL*(CPS*TMP*TREF/DES+HLAM)/ (CPL*TREF)+ 
1TMP 

IF(H(1,M) .GE.HLMP) GO TO 630 
IP=3 

GO TO 623 

630 DO 650 1=2, NOW 

IFCICOUNT.NE . 1) GO TO 635 

CALL CHANGE (HO(I-l,M) ,HO(I ,M) ,H0(I+1 ,M) ,HO(I ,M) ,HO(I ,M-1) ,IPO, 

1 AKO 1 , AK02 , AK03 , AK04 ) 

C0NST(I,M)=~AK01*X1*(T0(I-1,M)-T0(I,M))+AK02*X1*(T0(I,M)-T0(I+1,M) 
l)+2 .*AK04*(T0(I ,M) -TO(I ,M-1))/YL1 
635 CALL CHANGE (H( I -1,M) ,H(I ,M) ,H(I+1 ,M) ,H(I ,M) ,H(I ,M-1) ,IP,AK1 , 
1AK2.AK3.AK4) 

VK 1 = ( AK 1+AK2 ) *X 1 
VK2=(AK3+AK4)/YLi 
VK3=2 .*AK3/YL1 
VK4=2.*AK4/YL1 
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GO TO(643 , 645 , 647 ) , IP 
C ICE 

643 H(I»M)*1./(1. +VAS* ( VK1+VK4+VA9 ) )*(H0 ( I , M) -VA3* ( -XI* (AK2*T( 1+1 , M)+ 
1AK1*T(I-1,M) ) -VK4*T(I,M-1)+C0NST(I ,M) -VA9*(G2TIME*2 . -T0(I ,M) ) ) ) 
T(I ,M)=H(I ,M)/ (TREF*CPS) 

IF(H(I,M).LE.HSMP) GO TO 650 
IP 31 1 

C MELTING POINT 

645 H(I,M)=HG(I ,M) -VA3*( (VK1+VK4+VA9)*TMP-X1*(AK2*T( 1+1 ,M)+AK1* 
1T(I-1,M))-VK4*T(I »M-1)+C0NST(I ,M)-VA9*(G2TIME*2. -TO(I ,M) )) 
T(I,M)=TMP 

IF(H(I,M) .GT.HSMP.AND.H(I,M) .LT.HLMP) GO TO 650 
IF (IP -EQ. 1 , OR. IP.EQ. 3) GO TO 650 
C LIQUID WATER 

647 H(I ,M)=1 . / ( 1 . +VAL* (VK1+VK4+VA9 ) )* (HO C I ,M)-VA3*( (VK1+VK4+VA9)*VA6- 
1X1* (AK2*T ( 1+ 1 , M ) +AK 1*T ( I - 1 , M ) ) - VK4*T ( I , M- 1 ) +CONST ( I , M) - VA9* 
2(G2TIME*2. -TO(I ,M)))) 

T(I ,M)=H(I ,M)/ (TREF*CPL) -DEL* ( CPS*TMP*TREF/DES+HLAM ) / (CPL*TREF ) + 
1TMP 

IF(H(I,M).GE.HLMP) GO TO 650 
IP=3 

GO TO 643 
650 CONTINUE 

IF(ICOUNT.NE.l) GO TO 655 

CALL CHANGE (HO (N- 1 ,M) ,HO(N,M) ,H0(N-1,M) ,HO(N,M) ,H0(N,M-1) , IPO, 
1AK01 ,AK02 , AK03 , AK04) 

CONST (N ,M) =-AK01*2 . *X1* (TO(N- 1 , M) -TO (N , M) )+2 . *AK04* (TO (N . M) - 
1T0(N,M-1))/YL1 

655 CALL CHANGE (H(N-1,M) ,H(N,M) ,H(N-1,M) ,H(N,M) ,H(N,M-1) ,IP,AK1,AK2, 
1AK3,AK4) 

VK1=(AK1+AK2)*X1 
VK2=(AK3+AK4)/YL1 
VK3=2.*AK3/YL1 
VK4=2 . *AK4/YL1 
GO T0(663,665,667),IP 
C ICE 

663 H (N , M)=l . / ( 1 . +VAS* (VK1+VK4+VA9 ) )* (HO (N , M) -VA3* ( -X1*AK1*2 .*T(N-1.M) 
1 - VK4*T(N , M- 1 ) +CONST (N , M) -VA9* (G2TIME*2 . -TO(N , M) ) ) ) 
T(N,M)=H(N,M)/(CPS*TREF) 

IF(H(N,M) . LE.HSMP) GO TO 670 
IP=1 

C MELTING POINT 

665 H(N,M)=HO(N,M) -VA3* ( (VK1+VK4+VA9 )*TMP-X1*AK1*2 . *T(N- 1 , M) -VK4* 

1T(N , M- 1 )+CONST (N ,M) - VA9* (G2TIME*2 . -TO (N ,M) ) ) 

T(N,M)=TMP 

IF(H(N,M) .GT.HSMP.AND.H(N,M) .LT.HLMP) GO TO 670 
IF(IP.EQ. 1 .OR. IP.EQ. 3) GO TO 670 
C LIQUID WATER 

667 H(N,M)=1 . / (1 .+VAL*(VK1+VK4+VA9) )*(HO(N,M) -VA3* ( (VK1+VK4+VA9 )*VA6- 
1X1*2 .*AK1*T(N-1 ,M) -VK4*T(N,M-1)+C0NST(N,M) -VA9*(G2TIME*2 . - 
2T0(N,M)))) 

T (N , M) =H(N , M) / (TREF*CPL) -DEL* (CPS*TMP*TREF/DES+HLAM) / ( CPL*TREF)+ 
1TMP 

IF(H(N,M) .GE.HLMP) GO TO 670 
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60 TO 663 
670 GO TO 677 
672 DUMAX=0 . 

CHECK THE CONVERGENCE FOR CONVECTION BOUNDARY CONDITIONS 
OF ONE LAYER CASE ONLY 

DO 675 1=1 ,N 
DO 675 J=2 , MOW 

DIF(I,J)=ABS(T(I,J)-PHOLD(I,J)) 

DUMAX=AMAX1(DIF(I,J), DUMAX) 

675 CONTINUE 
GO TO 6S0 

CHECK CONVERGENCE OF THE GAUSS -SEIDEL ITERATION 

677 DUMAX=0 . 

DO 678 1=1, N 
DO 678 J=1,MJ 

DIF(I,J)=ABS(T(I,J)-TOLD(I,J)) 

DUMAX=AMAX1 (DIF( I , J) .DUMAX) 

678 CONTINUE 

CHECK SEE IF NEED ANY MORE ITERATION 

680 IF (DUMAX. LE.EPSI) GO TO 682 
I C0UNT=I C0UNT+ 1 
IF(L.EQ.l) GO TO 260 

IF((IBC1.EQ. 1) .AND. (IBC2.EQ. 1))G0 TO 315 
IF((IBC1.EQ.2) .AND. (IBC2.EQ.2))G0 TO 300 


682 TITAU=TAU*3600 . 

DETERMINATION OF WHETHER ICE LAYER SHOULD BE SHED 

IF(ISH.EQ.l) GO TO 685 
IF(JSH.EQ.3) GO TO 685 

IF((JSH.EQ.2) .AND. (NIP.LT.NFSP)) GO TO 685 
HINC=H(1,INC) 

HINC=HAFP+(HINC-HSMP)/2 . 

IF (HINC . LT . 2IAFP) GO TO 685 
JSH=2 

IF(IRP.EQ.l) JSH=3 
NISP=NIP 


DETERMINATION OF WHETHER THE PROGRAM SHOULD BE TERMINATED 
BY PHYSICAL TIME 

685 TCOUNT=0 .015 

(TAU . GE . TCOUNT) GO TO 999 

DETERMINATION OF WHETHER THE PROGRAM SHOULD BE TERMINATED 
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C BY INPUT DATA OF TIME STEPS 
C * 

IF(NIP.GT.NTSP) GO TO 999 

DETERMINATION OF WHETHER TEMPERATURES AND ENTHALPYS 
FOR THE ICE LAYER SHOULD BE STORED FOR THE PURPOSE OF 
RECREATING ICE AND SHEDDING THE ICE 

IF(IRF.EQ.l) GO TO 694 
IF(NIP.LT.NFSP) GO TO 694 
IF(IRP.EQ.l) GO TO 694 
IFC(ISH.EQ.l) .AND. (JSH.EQ.l)) GO TO 694 
687 DO 692 1=1, N 
DO 692 J=INC,M 
T(I, J)=TIN 

IFCTIN.LT. TMP) GO TO 690 

H C I , J)=CPL*TREF*T(I , J ) +DEL* (CPS lfr TMP' ,f TREF/DES+HLAM ) -TMP*TREF* 
1CPL 

690 H(I,J)=CPS*T(I,J)*TREF 
692 CONTINUE 
IRP=1 

694 IF(ISH.EQ.l) GO TO 69S 
IF((JSH.EQ.2) .AND. (NIP.LT.NFSP))MJ=INC 
IF ( ( JSH . EQ . 2 ) . AND . (NIP . GE . NFSP) )M J=M 
IF ( ( JSH . EQ . 3 ) . AND . (NIP . GE . NFSP) )MJ=INC 
IF(JSH.EQ.l)MJs*J 

695 IF(LA.NE.IFREQ) GO TO 200 
TITAU=TAU*3600 . 

PRINT OUT OF PROGRAM 

WRITE (10, 766 )TITAU 
DO 697 1=1, N 
DO 697 J=1,MJ 
TR(I, J)=T(I, J) 

TE(I , J)=TR(I , J)*TREF 
697 CONTINUE 

WRITE(I0,768) ( (TE(I, J) , 1=1, N) , J=1,MJ) 

IF ( IRF . EQ . 2 . AND . NIP . EQ . NFSP ) GO TO .699 
WRITE (10,770) ICOUNT 

699 LA=0 
GO TO 200 

700 FORMAT (80A1) 

701 F0RMAT(56X,I3/56X, I3/S6X,I3/56X,F10.6/56X,F10.6/56X,F10.6) 

702 FORMAT(7X,I2,7X,F10.5,9X,F13.6,14X,F13.6) 

703 F0RMAT(53X,I3/53X,I3/53X,I3/53X,I3) 

704 FORMAT(53X, I3/53X, I3/53X, I3/53X, 13) 

705 F0RMAT(53X,F11.5/53X,F11.5/53X,F11.5/53X,F11.5) 

706 F0RMAT(53X,I3/53X,I3) 

707 FORMAT (46X ,F13 . 6/46X ,F15 . 6/46X ,F13 . 6/46X ,F15 .6) 

708 F0RMAT(46X,F13 ,6/46X,F13 .6) 

709 F0RMAT(61X,I3/47X,F13.6/47X,F13.6/47X,F13.6/47X,F13.6/47X,F13 .6 
1/47X,F13 . 6/47X.F13,. 6) 

710 F0RMAT(47X,F13.6) 
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711 FORMATC50X,F13.6/50X,I3/50X,F13.6/50X,I3/50X,F13.6) 

712 FORMAT(50X,I3) 

713 F0RMAT(68X,I2/68X,I2/58X,I3/50X, I3/50X,I3/50X,I3) 

716 F0RMAT(46X,F13.6/46X,F13.6) 

717 FORMAT (/ ' TOTAL NUMBER OF SLABS’ ,38X,'L=',I3/ ' TOTAL NUMBER 
10F NODES IN THE X-DIRECTION OF THE GRID' ,3X, 'N=' ,13/ ’ TOTAL 

2 NUMBER OF NODES IN THE Y-DIRECTION OF THE GRID* , IX, 'M=' , 

3 13/ ' TOTAL LENGTH OF COMPOSITE SLAB IN THE Y-DIRECTION 1 , 

4 6X, 'TLEN1=' ,F13.6, ' INCHS'/ ' TOTAL LENGTH OF THE HEATER IN THE 
5X-DIRECTION' ,8X,'TLEN2=',F13.6,’INCHS7 ' GAP WIDTH' ,46X, 

6'TLEN3=' ,F13.6,'INCHS') 

718 FORMATC INTERNAL HEAT GENERATION IN SLAB NUMBER', 9X,' IJ=' , 
1I3/33X, 'BETWEEN NODE N01=* .I3/36X, 'AND NODE N02=’,I3) 

719 FORMAT ( / / 10X , ' THERE IS NO HEAT SOURCE PRESENT ') 

720 FORMAT (/ ' POINT HEAT SOURCE IS PRESENT AT' , 17X, 'NODE=' ,13/ 

133X,' BETWEEN SLAB L1=\I3/36X, ' AND SLAB L2=’,I3) 

721 FORMAT (/ 1 CONSTANT HEAT INPUT OF',28X,'Q = \F13 .6, ’WATTS/IN*IN' ) 

722 FORMATC/ ' ON-TIME FOR HEAT INPUT' ,27X, 'TON=' ,F13 ,6, ' SECS'/ 

1' OFF-TIME FOR HEAT INPUT' ,25X, 'TOFF=' ,F1 3.6, ’ SECS') 

724 FORMATC/ ’ STEP FUNCTION HEAT INPUT OF' ,28X, 'QV=' ,F13.6, 

1' WATTS/ INKIN' ) 

726 FORMATC/ ' RAMP FUNCTION HEAT INPUT OF' ,28X, 'QV=' ,F13.6, 

1' WATTS/ IN* IN') 

727 FORMATC/ ' SINE FUNCTION HEAT INPUT OF’ ,28X, 'QV=' ,F13.6, 

1’ WATTS/ IN*IN') 

729 FORMATC// ' CONSTANT TEMPERATURE AT X=l',17X,’ TX1=' .F13.6, 

1’ DEG . F ' ) 

730 FORMATC// ' CONVECTION OCCURS AT X=1 ’/I IX, ’AMBIENT TEMPERAT 
1URE' ,5X, *TG1=' ,F13.6, 'DEG.F'/llX, 'HEAT TRANSFER COEFF. ' ,5X, 'HI =’ , 
2F15.6, 'B.T.U/HR.FT.FT.DEG.F') 

732 FORMATC// ' CONSTANT TEMPERATURE AT X=2' ,17X, 'TX2=* ,F13.6, 
l'DEG.F') 

734 FORMATC// ’ CONVECTION OCCURS AT X=2 '/I IX, 'AMBIENT TEMPERAT 

1URE' ,5X, 'TG2=' ,F13.6, 'DEG.F’/I’X, 'HEAT TRANSFER COEFF. ' ,5X, 'H2 = ' , 
2F1S. 6, 'B.T.U/HR.FT.FT.DEG.F' ) 

736 FORMATC/ ' THE INITIAL TEMPERATURE IN THE COMPOSITE SLAB TIN =' , 
1F13.6, 'DEG.F'// ' THE REFERENCE TEMPERATURE 20X , ’ TREF =',F13.6, 
2'DEG.F') 

738 FORMATC// ' THE PHASE CHANGE IN THE ICE LAYER IS CONSIDERED ') 

740 FORMATC/ ' LATENT HEAT OF ICE ',21X,' HLAM =’ ,F13.6, 'B.T.U./LB'/ 

1' THERMAL CONDUCTIVITY OF WATER' , 12X, 'AKL = ' ,F13.6, 'B.T.U./HR. 
2FT. DEG.F'/ ’ THERMAL DIFFUSIVITY OF WATER ',12X,'ALPL =',F13.6, 

3' FT. FT/HR. ’/ ’ DENSITY OF WATER ',24X,'DEL =\F13.6, 'LB/CU.FT'/ 
4 ' SPECIFIC HEAT * DENSITY OF WATER \8X,’CPL =' ,F13.6, 'B.T.U./ 
5CU. FT. DEG.F'/ ' DENSITY OF ICE ',26X,'DES = \F13. 6, 'LB/CU.FT. ’ ) 
750 FORMATC/ ' SPECIFIC HEAT * DENSITY OF ICE '.lOX.'CPS =',F13.6, 

1’ B.T.U. /CU. FT. DEG.F') 

752 FORMATC/ ' ICE MELTING TEMPERATURE' ,18X,'TMP =' ,F13.6, 'DEG.F' ) 

753 FORMATC// ' THE PHASE CHANGE IN THE ICE LAYER IS NOT CONSIDERED’) 

754 FORMATC/ ' INITIAL TIME STEP' ,26X, 'DTAUI =' ,F13.6, 'SECS'/28X, 
l'FOR TIME STEPS NISP =',13/ ' INTERMIDIATE TIME STEP',21X, 
2’DTAUM =' ,F13 . 6 , ' SECS ' /28X, ' FOR TIME STEPS NMSP =',13/ 

3' FINAL TIME STEP’ ,28X, 'DTAUF =' ,F13.6, ’SECS') 

756 FORMATC/// ' SHEDDING ICE LAYER IS CONSIDERED') 
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758 FORMAT(/ ' INITIAL TEMPERATURE IS CHANGED FOR TIME STEPS NFSP = ', 
113/ 1 THE CYCLE FOR SHEDDING' ,24X, ' ICL =’,13/ ' FOR SHEDDING 
2 AT 2ND CYCLE FOR TIME STEPS NHSP-' ,13) 

760 FORMAT (' FREQUENCY OF TIME STEP/PRINT OF OUTPUT', 6X, 

l'IFREQ=' ,13/ ' FOR TIME STEPS STOP THE PROGRAM' , 13X, 'NTSP =',I3) 
762 FORMAT(///28X, ' TEMPERATURE PROFILE IN DEGREES F ') 

764 F0RMAT(/2X, 'HEATER START AT NODE IN X-DIRECTION' ,SX, 'NODEG*' ,13) 
766 F0RMAT(///35X, ' TIME TAU =’,F1S.6,’ SECS') 

768 F0RMAT(/5X,20F6.2) 

770 F0RMATC/2X, 'PASS®' ,13) 

775 FORMATC ' ,80A1) 

999 STOP 
END 


SUBROUTINE FOR CONSTANT STEP FUNCTION HEAT SOURCE 


FUNCTION QHEAT2 (TAU , DTAU, TON , TOFF ,QV) 
TAUR=TAU 

TAUR=TAUR+0 . 5*DTAU 
IN=IFIX(TAUR/ (TON+TOFF) ) 

IP=IN+1 

B =IN*TOFF+IP*TON 
C =IP* (TON+TOFF) 

QHEAT2=QV 

IF( (B . LT . TAU) . AND . (TAU . LT . C ) )QHEAT2=0 . 

RETURN 

END 


SUBROUTINE FOR RAMP FUNCTION HEAT SOURCE 


FUNCTION QHEAT3 (TAU , DTAU , TON , TOFF , QV) 
TAUR=TAU 

TAUR=TAUR+0.5*DTAU 
IN=IFIX (TAUR/ (TON+TOFF) ) 

YP=TAUR- IN* (TON+TOFF) 

A=QV/TON 

IF(YP.EQ.O) GO TO 20 
IF (YP -TON) 10, 10, 20 
10 QHEAT3=YP*A 
GO TO 50 
20 QHEAT3=0. 

50 RETURN 
END 


SUBROUTINE FOR SINE FUNCTION HEAT SOURCE 


FUNCTION QHEAT4 (TAU , DTAU , TON , TOFF , QV) 
TAUR=TAU 
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YP=TAUR- IN* (TON+TOFF) 

XP=YP*2.*3 . 14159/ (TON+TOFF) 

QHEAT4=5 .*SIN(XP)+QV 

RETURN 

END 


SUBROUTINE FOR CORRECTING THERMAL CONDUCTIVITIES FOR PHASE CHANGE 


SUBROUTINE CHANGE (HL1 ,HL2 ,HL3 ,HL4 ,HL5 , IPP,AK1 ,AK2 ,AK3 ,AK4) 
COMMON / SEA 1 /HSMP , HLMP » AKS , AKL 
DIMENSION HL(5) ,IP(5) ,CK(5) 

HL(1)=HL1 
HL(2)=HL2 
HL(3)=HL3 
HL(4)=HL4 
HL(5)=HL5 
DO 50 1=1,5 

IF(HL(I).LE.HSMP) GO TO 10 
IF(HL( I). GT.HSMP.AND.HL(I).LT. HLMP) GO TO^O 
IF(HL(I) .GE.HLMP) GO TO 30 
10 IP(I)°1 
CK(I)=AKS 
GO TO 50 
20 IP(I)«=2 

X=(HLMP-HL(I))/(HLMP-HSMP) 

CK(I)«X*AKS+(1 . -X)*AKL 
GO TO 50 
30 IP(I)=3 
CK(I)=AKL 
50 CONTINUE 
IPP=IP(2) 

AKl=(CK(l)+CK(2))/2. 

AK2=(CK(3)+CK(2))/2. 

AK3=(CK(4)+CK(2))/2. 

AK4=(CK(5)+CK(2))/2. 

RETURN 
END 

/* 

//G0.FT06F001 DD SYSOUT=A,OUTLIM=9990 
//GO.SYSIN DD * 


DATA FOR THE PROBLEM 


COMMENTS: 

IF PHASE CHANGE IS CONSIDERED IG=2 
IF PHASE CHANGE IS NOT CONSIDERED IG=1 

FOR CONSTANT TEMPERATURE BOUNDARY CONDITIONS AT X=1 IBCI=1 
" » .. .. " AT X=2 IBC2=1 

" CONVECTION BOUNDARY CONDITIONS AT X=1 IBC1=2 

" ” " " AT X=2 IBC2=2 
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IF NO HEAT SOURCE IS PRESENT IH * 1 
IF POINT HEAT SOURCE IS PRESENT IH * 2 
IF INTERNAL HEAT GENERATION IN A SLAB OCCURS IH >3 
IF CONSTANT HEAT INPUT IS USED I HQ * 1 
IF STEP FUNCTION HEAT INPUT IS USED IHQ = 2 
IF RAMP FUNCTION HEAT INPUT IS USED IHQ « 3 
IF SINE FUNCTION HEAT INPUT IS USED IHQ = A 


DATA FUR THE COMPOSITE BODY 

Vrrtrt VrVWnYiWWMr iWriVrtrttWrVrVrVf ifrrt 

TOTAL NUMBER OF SLABS U= 6 

TOTAL NUMBER OF NODES IN THE X-DIRECTION OF THE GRID N« 20 

TOTAL NUMBER OF NODES IN THE Y-DIRECTION OF THE GRID M= 29 

TOTAL LENGTH OF COMPOSITE SLAB IN THE Y-DIRECTION TLEN1* 0.413000 INCHS 

TOTAL LENGTH OF THE HEATER OF THE GRID TLEN2= 0.255000 INCHS 

GAP WIDTH TLEN3= 0.015000 INCHS 

DATA FOR EACH SLAB: 


NUMBER OF NODES 

LENGTH 

THERMAL CONDUCTIVITY 

THERMAL DIFFUSIVITY 


IN INCHS. 

IN B.T.U./HR. .FT. 'F. 

IN FT*FT/HR. 

S 

0.08700 

66.500000 

1.650000 

5 

0.05000 

0.220000 

0.008700 

2 

0 . 00400 

7.600000 

0.136000 

3 

0.01000 

0.220000 

0.008700 

3 

0.01200 

8.700000 

0.150000 

16 

0.25000 

1.416000 

0.049200 


DATA FOR THE HEATER: 


TYPE OF HEAT SOURCE 

FOR IH“2 POINT HEAT SOURCE AT 

BETWEEN SLAB 
AND SLAB 

FOR IH=3 INTERNAL HEAT GENERATION IN SLAB 

BETWEEN NODE 
AND NODE 

TYPE OF HEAT SOURCE USED 
CONSTANT HEAT INPUT OF 
ON-TIME FOR THE VARIABLE HEAT INPUT 
OFF-TIME FOR THE VARIABLE HEAT INPUT 
VARIABLE HEAT INPUT OF 

BOUNDARY AND INITIAL CONDITIONS i’ 


IH= 3 
NODE= 9 
LI- 2 
L2= 3 

IJ= 3 
N01= 9 

N02= 10 
IHQ= 2 

Q = 25 . OOOOOWATTS/IN*IN 

TON = 10 . OOOOOSEC . 

TOFF = 10.0000 SECS. 

QV = 25.00000WATTS/IN*IN 


TYPE OF BOUNDARY CONDITION 

IF IBC1=1, CONSTANT TEMPERATURE 
IF IBC2-1, " " 

IF IBC1=2, AMBIENT " 

HEAT TRANSFER COEFF 
B.T.U./HR.FT*FT*'F 
IF IBC2=2, AMBIENT TEMPERATURE 
HEAT TRANSFER COEFF. 




AT X=1 

IBC1= 2 



AT X=2 

IBC2= 2 

AT 

X=1 

TX1 = 

20.000000 DEG.F 

AT 

X«2 

TX2 = 

20.000000 DEG.F 

AT 

X=1 

TGI = 

-4.00000 DEG.F 

AT 

X=1 

HI = 

1.000000 

AT 

X=2 

TG2 = 

-4.00000 DEG.F 

AT 

X-2 

H2 = 

1000000.000000 


B.T.U/HR.FT.*FT*'F 

THE INITIAL TEMPERATURE IN THE SLAB 

THE REFERENCE TEMPERATURE 


TIN' = -4.000000 DEG.F 

TREF = 32.000000 DEG.F 
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PROPERTIES OF WATER: 


IS THE PHASE CHANGE CONSIDERED IF YES 
LATENT HEAT OF ICE 
THERMAL DIFFUSIVITY OF WATER 
" CONDUCTIVITY " " 

DENSITY OF WATER 

SPECIFIC HEAT^DENSITY OF WATER 

DENSITY OF ICE 

SPECIFIC HEAT*DENSITY OF ICE 

ICE MELTING TEMPERATURE 


IG=2 IF 

NO 

IG=1 IG= 2 

HLAM 

s 

143.40000 

B.T.U./LB 

ALPL 

s 

0,0051 

FT*FT/HR. 

AKL 

s 

0.32 

B.T.U./HR.FT.'F. 

DEL 

3 

62.40000 

LB./CU.FT 

CPL 

= 

62.21280 

B.T.U./CU.FT.’F. 

DES 

s 

57.40000 

LB./CU.FT 

CPS 

s 

28.81480 

B.T.U./CU.FT.'F. 

TMP 

B 

32.00 

DEG. 'F. 


.ADDITIONAL DATA: 


INITIAL TIME STEP DTAUI * 

FOR TIME STEPS NISP =75 

INTERMIDIATE TIME STEP DTAUM = 

FOR TIME STEPS NMSP =168 

FINAL TIME STEP DTAUF = 


O.IOOOOO SECS. 


0.001000 SECS. 


0.100000 SECS. 


FREQUENCY OF TIME STEP PF.R PRINT OF OUTPUT I FREQ = 3 
IS THE SHEDDING ICE LAYER CONSIDERED IF YES ISH=2 IF NO ISH=1 
DOES THE INITIAL TEMPERATURE CHANGE IF YES IRF=2 IF NO IRF=1 
INITIAL TEMPERATURE IS CHANGED FOR EACH TIME STEPS NFSP =300 
THE CYCLE FOR SHEDDING ICL = 2 

FOR SHEDDING AT 2ND CYCLE FOR TIME STEPS NHSP =463 
FOR TIME STEPS STOP THE PROGRAM NTSP =600 

/* 

// 


ISH= 

IRF= 
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TABLE X. THERMAL PROPERTIES OF SELECTED MATERIALS 


Material 

Thermal 

Conductivity 

Spec. 

Heat 

Dif fuslvity 

Density 


Btu 

W 

Btu 

J 

ft 2 in 5 m^ 

lb 

K& 

f£ 

-hr-°F a 

i~°C 

lb-°F 

Kg-°C 

hr 

■“ s 

ft 3 

m 3 

Aluminum 

(soft) 

Aluminum 

Alloy 

124 

215 

0.23 

963 

3.20 

8.26 

169 

2705 

75S-T6 

66.5 

115 

0.23 

963 

1.65 

4.26 

175 

2800 

Beryllium 

Copper 

60.5 

105 

0.10 

419 

1.17 

3.02 

515 

8250 

Copper 

220 

381 

0.092 

385 

4.29 

11.1 

557 

8900 

Nickel 

36 

62.3 

0.105 

440 

0.62 

1.59 

556 

8900 

Nichrome 

80/20 

7.6 

13.2 

0.107 

448 

0.138 

0.356 

515 

8250 

Stainless 

Steel 304 

Glass 

("E" 

8.7 

.15.1 

0.118 

494 

0.15 

0.387 

495 

7930 

Glass) 

0.6. 

1.0 

0.19 

796 

0.020 

0.052 

159 

2550 

Epoxy 

Resin 

Epoxy/ 

Glass 

0.10 

0.173 

0.25 

1050 

0.0053 

0.0137 

75 

1200 

filled 

32% 

0.22 

0.38 

0.23 

963 

0.0087 

0.0225 

110 

1760 

Neoprene 

0.108 

0.187 

0.295 

1235 

0.0041 

0.0105 

90 

1440 

Nylon 0.14 

Polytetra- 

fluoroethylene 

0.24 

0.4 

1675 

0.0050 

0.013 

70 

1120 

(Teflon) 

0.15 

0.26 

0.25 

1047 

0.0043 0.011 

138 

2210 
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TABLE 1. (CONTINUED) 


Material 

Thermal 

Conductivity 

Spec. 

Heat 

Diffusivity 

Density 


Btu 

W 

Btu 

J 

ftl ! 

0 5 2_ 

lb 

& 


ft-hr-°F 

m-°C 

lb~°F 

Kg~°C 

hr 1 

0 S 

ft 3 

m 3 

Poly- 

trifluoro- 

chloro- 

ethylene 

(KEL-F) 0.04 

0.07 

0.216 

904 

0.0014 0.0036 

130 

2080 

Water 

32°F 

0.320 

0.554 

0.997 

4174 

0.0051 0.0132 

62.4 

999.6 

Ice 
(Pure) 
32 °F 

1.293 

2.238 

0.5057 

2117 

0.0445 0.115 

57.2 

916.3 

14°F 

1.356 

2.347 

0.5038 

2109 

0.0469 

0.121 

57.3 

917.9 

-4°F 

1.416 

2.451 

0.5020 

2102 

0.0492 

0.127 

57.4 

919.5 


Latent Heat of Fusion = 143.4 Btu/lb (333.6 KJ/Kg) 


TA3LE 2. HEATING TIMES REQUIRED TO PRODUCE GIVEN TEMPERATURE RISES 
OVER THE CENTERLINES OF THE HEATER ELEMENT AND THE GAP 
BETWEEN ELEMENTS FOR FOUR GAP WIDTHS. 


fa 

(0 

u 


00 


• • 
lO lO 


00 

* 

fa 


o\ 


© tf) CM 

• . • 

H f-4 xT 


fa 

o 

vO 


t 4 






0 ) 

00 

*H 

o 

CO 

H 

u 

• 

• 

• 

• 

« 

(0 

a) 

ES 

m 

vO 

V 0 

VO 

vO 


H 


m 

* 

VO 


fa 

o 

vO 


CM 


a 

o 

•H 

AJ 

u 

8 

44 

<0 

8 

U 






M 


«•*% 






/*N 

|| 





a> 

in 

m 

in 

in 

m 

in 

m 

m 






44 

• 

• 

» 

• 

« 

9 

9 

9 






CO 

CM 

CM 

CM 

PM 

CM 

CM 

CM 

CM 



a) at a> 



CD 


w 


w 




w 

•C 

i-i 

u A « H 



33 










ai 

co *h m ai 













<u 

u 

a s as 












/■s 

Vi 

e s s c /5 














80 S 













(0 

M H M CO 












w 

CO 

in CD 













0 ) 

10 • 00 0 ) 













fH 

co 0 co *3 













a 

<3 co e 





/-s 

t 



/•“S 




•3 

»H »rl fl 



fa 

CM 

CM 

in 

in 

OV 

Ov 

m 

in 


CO 

tf CO 0 CO 



<0 

9 

• 

• 

9 

9 

9 

9 

9 


u 

fa. CO fa. 44 



tf 

fH 

H 

<h 

aH 


H 

<M 

CM 


m 










w 




<r 

0 .§ 0 <r 

■3 











O 

o 

fa O faO 

fa 

fa 










O 

m 

u *rt b) cn 

CO 

0 










CM 


Jo 

iJ 

00 











2 

S H S S 5 

14 

* 










|J 

o 

0 000 

CO 

0 











p 

CM O H h m 

S 

H 

*4 









CM 

H 

O *4 O 0 H 



<0 






/**> 




. 

• (!)••• 

m 


44 

O . 

O 

0 

O 

O 

O 

O 

O 


O 

O N O O O 

CM 


CO 

• 

9 

• 

• 

9 

9 

9 

9 






CD 

H 

rH 

tH 

H 

*H 

H 

H 

fH 

A 





, 33 


'w' 




w 


w 

O 


u 

at 

u 


l 

a) 

a 


a a 

o o *o 

44 *H H 

m quill 

(0 3 (0 *H 

H I H ,fi 

3 3 3 W 

0 ) in h (j) >4 

u a u a a <u 

H W O fa 

>4 "H <0 
CD H 0) *4 

U Cl O 

III U M Cl 


fa 

+1 


2 

U |4 

to a) 

43 - 

3 
v> 


a 

i 

o 


cu 

V) 


II 

•H 


M 

a) 


43 
fa *4 


a ns ^ 


OH 
13! W 


o 

H 

o 

o 


o 

n 

o 


o 

in 

o 


o 

fa 

o 


97 




Cl 


n 


fa 

00 

00 

vO 

vO 

m 

4 -^ 

vD 

VO 

00 

I) 

fj 

co 

9 

9 

9 

9 

• 

# 

9 

• 

CM 


O 

CM 

CM 

n 

n 

•» 

<r 

m 

in 

w 

V 5 

q 

1 * 


L. 


I 


0 . ■■***>-■ ■ T 


TABLE 3. EFFECT OF VARIOUS ABRASION SHIELD MATERIALS 
ON ABRASION SHIELD-ICE INTERFACE TEMPERATURE 


Case 

Abrasion 

Shield 

Time 

Over 

to Reach 32 °F 
Center of Heater 
(Sec.) 

Time 

Over 

to Reach 32 °F 
Center of Gap 
(Sec.) 



[5] 

Present 

[5J 

Present 

l 

0.010" S.S. 

0.81 

0.9 

3.73 

4.2 

2 

0.010" 

nickel 

0.80 

0.9 

3.37 

3.8 

3 

0.010" 
nickel over 
0.010" S.S. 

1.10 

1.2 

3.36 

3.7 

4 

0.020" 
nickel over 
0.010" S.S. 

1.35 

1.5 

1.97 

2.2 

5 

0.010" S.S. 
over 0.010" 
nickel 

1.12 

1.2 

1.59 

1.9 

6 

0.010" S.S. 
over 0.010" 
Epoxy/Glass 
over 0.010" 
nickel 

2.12 

2.2 

2.34 

2.6 


Other Data as Follows: 


Substrate: 

Inner Insulation: 
Heater Element: 

Other Insulation: 

Gap between Elements: 


0.100" thick, 304 Stainless Steel 
0.020" thick, Epoxy/Glass Laminate 
Zero thickness, 0.510" wide 
0.010" thick, Epoxy/Glass Laminate 
0.050" wide 


98 



OKJGiNAi, 

OF POOR 


PAG* ($ 

QUALITY 


h b2’ T b2 



AIR 


ICE 


ABRASION SHIELD 

OUTER INSULATION 
HEATER ELEMENT 
INNER INSULATION 

SUBSTRATE 


99 


HEATER 


OWflWAL wwu MJ 

OF POOR QUALITY 


as 

H 

< 


U 

<; 




FIGURE 2. TWO-DIMENSIONAL FINITE DIFFERENCE REPRESENTATION 


FIGURE 3 
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OF POOR QUALITY 



Interior Mesh Point at the 
Center of the Gap 
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|»ix 1 



Interior Mesh Point at the 
Center of the Heater 



Interfacial Point 


FINITE DIFFERENCE GRID AT SELECTED POINTS 
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a. Substrate - Ambient Interface 



ambient 
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b. Abrasion Shield - Ambient Interface 
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c. Ice-Ambient Interface 


FIGURE 4. FINITE DIFFERENCE GRID AT 
SELECTED POINTS 
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thermal conductivity 
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FIGURE 7. EFFECT OF POWER DENSITY ON DE-ICER PERFORMANCE 



25 Watts/In 



FIGURE 8. COMPARISON OF RESPONSE OF POINT AND FINITE 
THICKNESS HEATERS 




Density 25 Watts/ in 

ate 0.087” 75S-T6 Alum 

Insulation ”N" Epoxy/Glass 

0.004" Nichrome 



FIGURE 9. COMPARISON OF DE-ICING TIME FOR VARIOUS RATIOS OF 
INNER TO OUTER INSULATION THICKNESS 
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FIGURE 10. EFFECT OF INSULATION THICKNESS RATIO AND THERMAL 
MATERIAL ON DE-ICING TIME 
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FIGURE 12a. COMPARISON OF INTERFACIAL TEMPERATURE RESPONSE FOR RAMP 
FUNCTION AND CONSTANT POINT HEAT SOURCE 
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FIGURE 13. EFFECT OF INITIAL ICE LAYER THICKNESS ON DE-ICER PERFORMANCE 
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FIGURE 14. TIME DIFFERENCE BETWEEN INTERFACE ABOVE CENTER OF HEATER AND 

INTERFACE ABOVE CENTER OF GAP TO ACHIEVE A GIVEN TEMPERATURE RISE 



De-Icer Pad Construction is Same 

Curve As Presented in Figure 12 
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FIGURE 15. TEMPERATURE PROFILES FOR THE TWO-DIMENSIONAL DE-ICER PAD 



De-Icer Pad Construction 

is Same as Presented in Figure 14 

Gap Width 0.050!' 
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FIGURE 16. EFFECT OF ICE MELTING ON INTERFACE TEMPERATURE 
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FIGURE 17. TEMPERATURE PROFILES WHEN ICE LAYER IS SHED, 1st CYCLE 
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FIGURE 18. TEMPERATURE PROFILES WHEN ICE LAYER IS SHED, 2nd CYCLE 







FIGURE 19. TEMPERATURE PROFILES WITH REFREEZING 





